If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Department of Mathematics, University of California, Riverside, California
Interdisciplinary Center for Quantitative Modeling in Biology, University of California, Riverside, California
Hair follicles (HFs) are mini skin organs that undergo cyclic growth. Various signals regulate HF cell fate decisions jointly. Recent experimental results suggest that transforming growth factor beta (TGF-
) exhibits a dual role in HF cell fate regulation that can be either anti- or pro-apoptosis. To understand the underlying mechanisms of HF cell fate control, we develop a novel probabilistic Boolean network (pBN) model on the HF epithelial cell gene regulation dynamics. First, the model is derived from literature, then refined using single-cell RNA sequencing data. Using the model, we both explore the mechanisms underlying HF cell fate decisions and make predictions that could potentially guide future experiments: 1) we propose that a threshold-like switch in the TGF-
strength may necessitate the dual roles of TGF-
in either activating apoptosis or cell proliferation, in cooperation with bone morphogenetic protein (BMP) and tumor necrosis factor (TNF) and at different stages of a follicle growth cycle; 2) our model shows concordance with the high-activator-low-inhibitor theory of anagen initiation; 3) we predict that TNF may be more effective in catagen initiation than TGF-
, and they may cooperate in a two-step fashion; 4) finally, predictions of gene knockout and overexpression reveal the roles in HF cell fate regulations of each gene. Attractor and motif analysis from the associated Boolean networks reveal the relations between the topological structure of the gene regulation network and the cell fate regulation mechanism. A discrete spatial model equipped with the pBN illustrates how TGF-
and TNF cooperate in initiating and driving the apoptosis wave during catagen.
Significance
Hair follicles (HFs) have become an emerging model system in the research of stem cell and wound-induced regeneration. We developed a literature-derived, single-cell-data-refined model on HF gene regulation dynamics and used it to investigate questions in HF cell fate regulations. We propose a new approach of synthesizing modeling and biological data that can be extended and applied to other systems in cell and developmental biology. Our model investigates the signaling regulation mechanisms underlying HF cell fate decisions, which provides new insights in the research of HF biology and cell fate regulation.
Introduction
Hair follicles (HFs) are stem cell-rich mini-organs in skin that can undergo oscillation-like cycles or regeneration throughout their lifetimes. Morphologically, the HF growth cycle includes three consecutive phases: anagen, the growing phase; catagen, the degenerating phase; telogen, the resting phase. Upon the telogen-to-anagen transition, HF stem cells are activated, which leads to the downward expansion of the HF until it reaches its maximum length. The HF will stay at this maximum length during the anagen while keeping active cell proliferation, maintaining the HF functions and producing the hair shaft. The most dynamic part of an anagen HF is at the bottom, usually referred to as the hair bulb. The center of the hair bulb consists of a group of specialized fibroblasts: the dermal papillae (DP), serving as the signaling headquarter (Fig. 1A). The DP is enveloped by a group of epithelial transient amplifying cells, the matrix (Mx) cells, which are fast dividing during the anagen phase. As Mx cells are pushed upward by cell proliferation, they differentiate into different components. Three major concentric layers are formed: medulla (MED), cortex/cuticle (CX), and inner root sheath (IRS), from inside to outside. Enclosing the IRS is the outer root sheath (ORS), which is directly derived from the HF stem cells located at the top of the HF, and the DP is connected to the connective tissue sheath that surrounds and protects the HF epithelium. Anagen may span days in mice, as opposed to humans, in which it may span a period of months or even years. At the end of anagen, an apoptosis wave is initiated from the bottom of the HF and propagates upward, leading to an upward involution of the HF epithelium. This transient degenerating phase is catagen. When the apoptosis wave stops, leaving the HF in its minimal form, the HF returns to the quiescent telogen phase. HF stem cells survive the coordinated apoptosis and await to be activated again upon the initiation of the next anagen.
Figure 1HF and its epithelial GRN. (A) An illustration of the bottom part of an HF. (B) The HF epithelial GRN summarized from literature and refined from scRNA-seq data, with scores from the PB method. To see this figure in color, go online.
In recent years, HFs have emerged as a leading model system for studying general mechanisms of stem cell regulation, tissue patterning, wound-induced regeneration, and tissue aging. Experimental results have uncovered how various signaling pathways cooperatively regulate cell division, differentiation, and apoptotic cell death in different compartments of the HF, which is crucial in regulating the HF cyclic growth and maintaining the HF functions (
). In particular, recent experimental studies have implied the seemingly dual roles—anti-apoptotic versus pro-apoptotic—of transforming growth factor beta (TGF-
) in regulating HF epithelial cell fate: upon the telogen-to-anagen transition, TGF-
2 counterbalances the refractory effects of bone morphogenetic protein (BMP) and activates anagen by activating the SMAD2/3 pathway (
); on the other hand, at the late stage of anagen, TGF-
2 activates the MAPK pathway and participates in inducing HF epithelial cell apoptosis, which terminates anagen and induces catagen (
). However, the regulatory mechanisms underlying the dual roles of TGF-
, as well as how TGF-
cooperates with other signals, are still unclear.
Although traditional reductionist approaches are powerful, they are limited by technological barriers and high experiment costs. As a result, combining modeling with empirical approaches to study tissue development and regeneration is useful and has been gaining interest. While many models have been developed and analyzed for HF morphogenesis during embryonic development (
). These HF cyclic growth models provide great insights on the interactions between signaling dynamics and HF development and growth control, yet the signaling events are either modeled using stochastic rules or an activator-inhibitor modeling framework, which cannot reflect the intra-cellular signaling transduction dynamics that underlies the HF cell fate decision mechanism. On the other hand, there has been an increase in modeling research on the signaling dynamics of apoptosis and TGF-
pathway. Many have chosen to adopt the Boolean network model, which is a popular modeling approach for studying signaling dynamics (
Network modeling of TGFβ signaling in hepatocellular carcinoma epithelial-to-mesenchymal transition reveals joint sonic hedgehog and Wnt pathway activation.
Boolean model of growth signaling, cell cycle and apoptosis predicts the molecular mechanism of aberrant cell cycle progression driven by hyperactive PI3K.
). However, to date, modeling studies particularly on the signaling regulation mechanisms in HF cells have been limited. Without a clear picture of the HF signaling transduction that is verified by experimental results, we cannot directly borrow the models originally developed from other systems and use them in HF modeling studies. In addition, many previous Boolean modeling studies focused more on the regulation mechanisms of TGF-
or apoptosis alone and rarely on the interactions of the TGF-
pathway with other signal pathways, despite that, at different stages of the HF growth cycle, various signals including TGF-
interact with each other to cooperatively regulate HF cell fate commitment.
In recent years, the fast-developed RNA sequencing techniques provide a much clearer picture of the skin and HF signaling events at the molecular level, and several datasets of skin and HF RNA sequencing data have become available (
). In parallel to the fast development of the RNA sequencing techniques, especially the single-cell RNA sequencing (scRNA-seq) technique, various computational biology methods that analyze these high-dimensional data have also emerged, which has strengthened the tie between biological data and mathematical models. Taking advantage of the scRNA-seq data and newly developed data-inference methods, in this paper, we first develop a Boolean network model from literature, then refine the model into a probabilistic Boolean network (pBN) model using recently available scRNA-seq data. Using this pBN model, we investigate the roles of TGF-
in regulating HF cell fate decisions in participation with BMP, tumor necrosis factor (TNF), and at different stages of the HF growth cycle. First, we propose a threshold-like switch in the TGF-
strength and demonstrate that this may necessitate the dual roles of TGF-
of anti- versus pro-apoptosis, as is suggested by experimental results. Next, we test our model to the high-activator-low-inhibitor theory of anagen initiation (
), with TGF-
serving as the activator and BMP the inhibitor. Our model results show concordance with the high-activator-low-inhibitor theory. Then, we use our model to test the roles of TGF-
and TNF in catagen initiation. This leads to our prediction that while both TGF-
and TNF can initiate HF epithelial cell apoptosis, TNF is more efficient than TGF-
. Additionally, we predict that shutting down pro-proliferation and pro-differentiation signals may further enhance apoptosis initiation. A discrete type of spatial model equipped with the pBN is also developed to test the two-step catagen initiation hypothesis and study the dynamics of apoptosis wave propagation during catagen. Finally, using our pBN model, we study the effects of gene knockout (KO)/overexpression (OE) on HF cell fate regulations, which might provide helpful predictions to guide future experiments. In addition, we apply the attractor analysis and stable motif analysis to the Boolean networks associated with the pBN model to help us understand how the topological structure of the network determines the nonlinear dynamics of the HF cell gene regulation.
Methods
The development of the HF pBN model mainly consists of two steps: 1) construct an initial Boolean network model for the HF gene regulation network (GRN) from literature, and 2) refine the Boolean network from published scRNA-seq data (
) and make it a pBN. We explain the details of the two steps next.
Constructing an initial HF GRN and the corresponding Boolean network from the literature
We first develop a GRN for HF matrix (Mx) cells, the group of transient amplifying epithelial cells at the bottom of an HF that envelop the DP (Fig. 1A). During anagen, Mx cells undergo fast division in response to various signals produced from DP, and on average it takes ∼12 h for an Mx cell to divide in mice, and ∼24 h in humans (
). Therefore, the HF Mx GRN, or equivalently the associated Boolean network, should have three major outcomes: cell division, differentiation, and apoptosis. For the signaling inputs of the GRN/Boolean network, at this point we include three signals: TGF-
, BMP, and TNF. TGF-
can either activate the cell division or apoptosis, BMP competes with TGF-
to initiate Mx cell differentiation, and TNF is known to play an important role in initiating cell apoptosis (
). The signaling transduction network is then built from the literature, including the KEGG PATHWAY Database:https://www.kegg.jp, shown in Fig. S1.
The initial GRN consists of two parts: a cell division/differentiation module and an apoptosis module (Fig. S1). In the cell division/differentiation module, TGF-
activates phospho-SMAD2/3, which then binds with SMAD4, activating further downstream signaling transduction that finally leads to cell cycling. BMP activates the cell differentiation signaling through the activation of phospho-SMAD1/5/8 signaling and its binding with SMAD4. For the cross-inhibition between TGF-
and BMP pathways, it is suggested that, in HF stem cells, p-SMAD2/3-SMAD4 target gene TMEFF1 negatively regulates BMP signaling (
). On the other hand, the BMP target genes inhibiting the TGF-
signaling in HF epithelial cells are not clear, yet cancer progression research suggests that CTGF and SMAD6 might be potential candidates: while SMAD6 may form the SMAD6/Smurf1 complex that inhibits the TGF-
signal, BMP may also inhibit CTGF, which is a potential activator of TGF-
signaling (
). However, we do point out that there are other research results showing that CTGF may antagonize BMP signaling, while SMAD6 can inhibit both TGF-
and BMP signaling (
). We will first include both CTGF and SMAD6 in the initial GRN, and later we will confirm their presence and roles from scRNA-seq data. Finally, SMAD7 may inhibit both TGF-
and BMP signaling (
Smad7 sensitizes tumor necrosis factor–induced apoptosis through the inhibition of antiapoptotic gene expression by suppressing activation of the nuclear factor-κB pathway.
), which is also included in the cell division/differentiation module. The apoptosis module is mostly developed based on both the KEGG PATHWAY Database and published GRN and Boolean modeling on apoptosis signaling (
). Nuclear factor κB (NF-κB) may be anti-apoptotic by activating BCL-XL/2 that inhibits mitochondria outer membrane permeabilization (MOMP), or by activating c-FLIP, which inhibits Caspase-8/10. In our Boolean network, to differentiate the signaling in the extrinsic versus intrinsic mitochondria pathway, we mathematically use two different nodes by placing or withholding an “m” in front of the gene’s name. For example, we have “BAX” for the extrinsic pathway and “mBAX” for the intrinsic mitochondria pathway. The intrinsic mitochondria pathway may lead to MOMP, which activates Caspase-9 and then Caspase-3/7. The two modules in the GRN are suggested to be linked through two pathways: ERK1/2 may inhibit SMAD4, and NF-κB may activate SMAD7, and both result in inhibitions of both BMP (SMAD1/5/8 and TGF-
) SMAD2/3 signaling and thus turn down cell division and differentiation. The activation of SMAD7 by NF-κB may in return inhibit the NF-κB activator IKK in the apoptosis module, thus acting as a self-inhibition on NF-κB.
With the initial GRN developed, we then develop its associated Boolean network model. Boolean network is a logical model that has been popularly adopted in systems biology, including the modeling study of signaling dynamics. In a Boolean network model, each node can take two possible values referring to the state of the node: 1 (ON) and 0 (OFF). The state of a node is updated by an associated Boolean function determined by the current states of its regulator nodes via the logic operators AND, OR, and NOT. We refer to (
) for a methodological introduction of the Boolean network modeling. To develop the Boolean functions associated with the initial HF GRN (Fig. S1), we assume that, when multiple regulators are involved, inhibitors always have a strong effect, thus inhibitors should always be included as AND NOT. On the other hand, we assume that activators are mostly supplement to each other, therefore we use OR unless an AND relation is clearly implied from biology. For example, in the model, MEK1/2 has two activators, RAF and TRADD, so we assume the following Boolean function:
while p-SMAD2/3-SMAD4 has two activators, p-SMAD2/3 and SMAD4, and clearly the Boolean function should be
Refine the Boolean network from scRNA-seq data into a pBN
) to refine our Boolean network model. First, we check if the signals in the initial GRN indeed express in the HF Mx cells. We find that PUMA, NOXA, and CTGF are not expressed; therefore, we remove these nodes and the edges linked to them from the model. We also find that SMAD8 and CASP10 are not expressed, yet in the model SMAD8 is grouped with SMAD1 and SMAD5, and CASP10 is grouped with CASP8, so the missing SMAD8 and CASP10 expressions do not affect our model.
Next, we use the pseudotime Boolean (PB) network inference method (
) to refine our Boolean model from the HF Mx scRNA-seq data. PB requires an initial Boolean network to start with, and it allows for the usage of the pseudotime order and binary expressions of cells to identify the most suitable Boolean functions, and the Boolean functions are scored based on how frequently they agreed with the pairs of input-output cells along the pseudotime trajectory. As shown in (
), Mx and its terminally differentiated progenies have four major subpopulations: germinative layer (GL), medulla (MED), cortex/cuticle (CX), and IRS. Diffusion map analysis clearly shows the three branches corresponding to MED, CX, and IRS that are derived from GL (Fig. S2, reproduced from data in (
)). We group the subclusters of the GL together with MED, CX, and IRS, depending on their adjacency to the three branches. Therefore, we obtain three trajectories where each of them starts from the GL and ends terminally differentiated, and we refer to the three trajectories as GL-MED, GL-CX, and GL-IRS, respectively (Fig. S2). We apply PB to each of the three trajectories with a few adjustments explained in Section S1. The scores along each trajectory as well as the average scores from the three trajectories are given in Table S1. Figs. S3–S5 show the expression heatmaps of the genes from our model along the three trajectories, reproduced from data in (
The ERK1/2 inhibition of SMAD4 receives straight 0 scores along all three trajectories, so we remove this edge from our Boolean network. This leaves the NF-κB→SMAD7→IKK pathway the only connection between the two modules in the GRN. In addition, after we remove the ERK1/2-SMAD4 edge, SMAD4 becomes another input node, and ERK1/2 becomes an end node that does not further affect the gene regulation dynamics; therefore, we remove the ERK1/2 node and the edges connecting to it as well. In our simulations, we will always keep SMAD4 ON except in the KO simulations. Finally, for each gene, we take the average score from the three trajectories, and, if multiple genes correspond to the same mathematical node, we take the highest average score among them (Table S1; Fig. 1B). We use these scores to further develop our model into a pBN model (
), where each Boolean function is associated with a probability given by the PB score. While a high PB score indicates confidence of the associated Boolean function, a low score implies the possibility of other regulations present in the system but absent from the current GRN. In this case, instead of searching for other potential regulators that would pose a great challenge due to various technical limitations, interpreting the score as a probability allows us to partially recover these uncertainty effects. Computationally, we use the general asynchronous method, where, at every computational step, 1) a node is randomly chosen to be updated, and 2) we randomly generate a random number
, compare it with the probability p of the Boolean function that updates this node, update this node if
, otherwise we skip it and choose the next random node to be updated. The computational scheme of updating the pBN model is shown in Fig. S6. We also model the signal inputs (BMP→SMAD1/5/8, TGF-
→SMAD2/3, RAS DAXX, and TNF→TRADD) in a similar probabilistic way to represent the strength of the signal inputs: in each computational step, we turn the signal input ON/OFF depending on its probability. In summary, the scRNA-seq data together with the PB inference method are used in three levels in the model refinement: 1) removing genes that are not expressed in HF Mx cells, 2) removing edges in the Boolean network that receive 0 scores, 3) using the PB scores for each Boolean function as the probability. The pBN modeling framework together with the general asynchronous method allows us to incorporate randomness into our model and recover the uncertainty effects.
A single run of the pBN model does not necessarily reach a steady state. In fact, even the associated deterministic Boolean network model shows oscillating dynamics in many cases, and cells can “commit” multiple fates at the same time (see Section S2); that is, more than one terminal node being ON at the same time (Fig. S7, slices with more than one fate). This implies that there might be other signals cross-inhibiting the pathways, yet they are not captured in the current GRN. To compensate for these missing factors in the systems, our model does not directly use the end gene nodes (p-SMAD1/5/8-SMAD4, p-SMAD2/3-SMAD4, CASP3/7) to indicate the cell fates. Instead, we mathematically add three modeling nodes (division, differentiation, and apoptosis) as terminal nodes representing each cell fate, and impose cross-inhibitions between them (dashed lined edges in Fig. 1B) to enforce the exclusive commitment of cell fates. With these cross-inhibitions among cell fates, cells can only turn ON one of the three cell states; yet we find another undetermined state (division = differentiation = apoptosis = OFF) (Fig. S9) in some simulations. Simulations and attractor analysis (see below for method details) show that, when all input signals are OFF, clearly this undetermined state is the only fixed point of the system; on the other hand, in the cases when the apoptosis module is triggered, this undetermined state results from the complex attractors in the system (Fig. S9). Mathematically, this undetermined state results from the AND NOT relations in the Boolean functions of the terminal fate nodes. At this point, we consider the undetermined state a modeled state due to possibly missing biological information in the current model, and it does not necessarily reflect a real biological state.
Now the development of our pBN is completed, altogether it has 38 nodes and 52 edges, and the Boolean functions are listed in Section S1. To model the average dynamics of the HF GRN from the pBN model, we take the average simulation results out of 10,000 simulations. Fig. 2 and Fig. S17 show an example result of the averaged dynamic simulations, and Fig. S18 shows the result from a typical single run of the pBN. We also refer to Section S3 for further computational assessment of the pBN and more discussions on the model setup.
Figure 2An example dynamic simulation from the pBN model, averaged from 10,000 simulations. Black = 1, White = 0. Initial conditions:
.
Attractor and stable motif analysis of the associated Boolean network
While the HF GRN dynamics is largely affected by the probabilities of the Boolean functions, the nonlinearity of the system dynamics also greatly depends on the graph topology of the Boolean network. Therefore, in addition to simulations of the pBN, we also apply the attractor and stable motif analysis (
) on the associated Boolean networks, which provide insights on the topological nature of the GRN dynamics. Details of the analysis methods and results are provided in Section S2 and Figs. S7–S12.
Results
A proposed HF cell fate regulation mechanism: TGF-
may regulate HF epithelial cell fates in a threshold-like switch fashion, jointly with BMP and TNF
The dual regulation effects of TGF-
have been observed in different systems. TGF-
2 can induce both the anti-apoptotic (SMAD2/3) or the pro-apoptotic (MAPK) pathway in airway epithelial cells, and studies have suggested this to be a threshold-dependent switch. When triggered by allergens, several airway cell types produce excessive TGF-
leading to increases in airway epithelial cells apoptosis (
). In immune systems, TGF-
triggers the differentiation of CD4+ cells to either induced Treg (iTreg) cells or T-helper q7 (Th17) cells, while iTreg cells produce signals to inhibit the differentiation to Th17 cells. Such a paradoxical regulatory role of TGF-
helps maintain the homeostasis of the system on the cell population level (
). In the HF system, experimental results of cultured HF stem cells suggest that TGF-
increases both the number and size of HF stem cell colonies below a certain concentration threshold (Fig. 4 from (
)). However, after passing the threshold, increasing TGF-
2 concentration will reduce the number and size of HF stem cell colonies. This seems to imply a threshold-like switch, with TGF-
2 activating the MAPK pathway when passing the threshold, similarly to that in the airway epithelial cells. Experimental results (Fig. 4 from (
)) also show that increasing the concentration of BMP4 decreases the number and size of HF stem cell colonies, while increasing both BMP4 and TGF-
2 does not change much of the number and size of HF stem cell colonies.
Simulation results from our pBN model with no TNF input (
) are shown in Fig. 3A–C, A′–C′. While BMP monotonically increases the chance of differentiation (Fig. 3B, B′) and monotonically decreases the chance of division (Fig. 3A, A′, arrow 2), TGF-
monotonically increases both the chances of cell division (Fig. 3A, A′, arrow 1) and apoptosis (Fig. 3C, C′), with the maximum chance of division ∼0.4. While BMP can initiate cell differentiation effectively, it is not hard to see that TGF-
is not very effective in initiating either division or apoptosis. We also observe that, when increasing both TGF-
and BMP along any contour line, the chance of cell division does not change (Fig. 3A, A′, arrow 3). Next, when a strong TNF input (
) is applied to the system, TNF significantly increases the chance of apoptosis despite the strengths of BMP and TGF-
(Fig. S19). With strong TNF (
), we find that TGF-
and BMP still monotonically increase the chances of division and differentiation, respectively. However, the maximum chances of division and differentiation have been greatly suppressed (<0.2) due to the elevated apoptosis.
Figure 3BMP and TGF-
regulation of cell fate. The y axis of the surface plots (A–F) and the colors in the heatmaps show the probability of that cell fate. (A, A′, B, B′, C, and C′) Cell fate (cell division, differentiation, or apoptosis) regulation by BMP and TGF-
in the pBN model with the uni-TGF-
node. DD′EE ′FF′) Cell fate regulation by BMP and TGF-
in the pBN model with separate TGF-
and Str-TGF-
nodes, the system dynamics shows a threshold switch. Simulations are taken with 0.1 increment on both TGF-
and BMP strength. Each group has 10,000 simulations and the results are taken at t = 1000. The x and y axis of each panel shows the strength of TGF-β or BMP, and the z axis of panels A– F shows the propability of the dedicated cell fate. To see this figure in color, go online.
The known dual roles of TGF-
in regulating HF epithelial cell fate—initiating HF stem cell cells division at the telogen-to-anagen transition, or initiating HF Mx cells apoptosis at the anagen-to-catagen transition—requires a more robust mechanism to differentially activate the SMAD2/3 or the MAPK pathways. The above-listed experimental evidence of a threshold-like switch in the TGF-
signaling transduction in both HF and airway epithelial cells suggests a potential regulatory mechanism. To install the threshold-based switch into the pBN model, we introduce another input node “strong TGF-
” (Str-TGF-β) into our pBN model. While both TGF-
and Str-TGF-
can initiate SMAD2/3, only Str-TGF-
initiates RAS and DAXX in the apoptosis module. See Section S1 for more details on the model setup.
With separate TGF-
and Str-TGF-
input nodes in the pBN model, simulation results show that, in the absence of TNF input, TGF-
first increases the chance of cell division, yet, after passing the threshold, it increases the chance of apoptosis while decreasing that of cell division (Fig. 3D, D′, arrow 1; Fig. 3F, F′), with a maximum chance of apoptosis ∼0.3 when the strongest Str-TGF-
input signal is applied. BMP still monotonically increases the chance of differentiation (Fig. 3E, E′) and monotonically decreases the chance of division (Fig. 3D, D′, arrow 2). In addition, when increasing both TGF-
and BMP along any contour line, the chance of division does not change (Fig. 3D, D′ arrow 3). All three findings (arrows 1, 2, 3 in Fig. 3D′) correspond to the experimental results of HF stem cell colony formation from (
). With the threshold switch of TGF-
, or, equivalently, separate TGF-
and Str-TGF-
input nodes in the pBN model, the introduction of a strong TNF input into the system causes the chance of apoptosis to be greatly elevated, accompanied with the suppression of both cell division and differentiation (Fig. S20).
Finally, using the attractor analysis on the corresponding Boolean networks, we investigate what leads to the multiple choice of cell fates in the long-term dynamic behavior of the system (Figs S9–S12). It turns out that whenever the apoptosis module is activated—either by TGF-
, Str-TGF-
, TNF, or their combinations—the system has a unique complex attractor resulting in multiple choices of cell fates. On the other hand, if only the cell division/differentiation module is activated, the system has only fixed points. Only in the case of TGF-
= BMP = ON and TNF = Str-TGF-
= OFF, the system has two fixed points, with one leading to cell division and the other differentiation. In all other cases with only the cell division/differentiation module activated, there is a unique fixed point in the system leading to a definite cell fate commitment: cell division or differentiation. Stable motif analysis shows that the only stable motifs are those corresponding to the fixed-point attractors, and there are no non-trivial stable motifs, defined as a motif that either includes only the input nodes and their state, or the whole fixed-point attractor state if one exists. We refer to Section S2 and Figs. S9–S12 for more details in the attractor and stable motif analysis.
Testing the model by an established theory: Anagen initiation requires both an increase in TGF-
and a decrease in BMP
Telogen, the resting phase of an HF, can be further divided into two sub-phases: the refractory and competent telogen (
). The refractory telogen is characterized by a high inhibitor/low activator profile, and it is followed by the competent telogen with a low inhibitor/low activator profile. The lowered inhibitor (represented by BMP) facilitates the anagen initiation when activator (including WNT and TGF-
) is elevated; on the other hand, with a high inhibitor level during the refractory telogen, even a high activator level cannot initiate anagen. We test these anagen initiation mechanisms using our pBN model with the threshold switch of TGF-
, that is, separate TGF-
and Str-TGF-
input nodes. When BMP is always high, the cell has little chance of division, and increasing TGF-
up to its threshold will not initiate division (division = apoptosis
, differentiation
; Fig. 4A). This means that, with high BMP during the refractory telogen, increasing TGF-
cannot initiate anagen. On the other hand, starting with a high BMP/low TGF-
profile, increasing TGF-
while decreasing BMP could increase the chance of division while decreasing the chance of differentiation (Fig. 4B), indicating the anagen initiation. In addition, no change in apoptosis is observed in this case (apoptosis
).
Figure 4Simulation results of the signaling regulations of anagen and catagen initiations. (A) During late telogen (BMP = ON, TNF = TGF-
= Str-TGF-
= OFF), increasing TGF-
alone cannot initiate anagen. (B) During late telogen, increasing TGF-
while decreasing Bmp may initiate anagen. (C–E) During late anagen (BMP = TGF-
= ON, Str-TGF-
= TNF = OFF), increasing Str-TGF-
alone (C), increasing TNF alone (D), or increasing both of them (E) may initiate catagen by increasing the chance of apoptosis. The y axis shows the probability of the dedicated signal or the cell fate. To see this figure in color, go online.
). We use our pBN model to examine the effects of TGF-
and TNF on apoptosis initiation. Starting with high BMP and high TGF-
but below the MAPK pathway initiation threshold, further increasing TGF-
alone increases the chance of apoptosis to ∼0.2 (Fig. 4C), while increasing TNF alone exhibits a more efficient apoptosis initiation with the chance of apoptosis >0.6 (Fig. 4D). Increasing both TNF and TGF-
beyond the threshold (Fig. 4E) shows a similar dynamic pattern with increasing TNF alone (Fig. 4D), implying that TNF dominates over TGF-
in the apoptosis initiation of a single cell. Our simulation results validate that, while both TGF-
and TNF alone can initiate apoptosis, TNF appears to be more efficient than TGF-
.
In addition to apoptosis initiation mechanisms, studies also suggest that, upon the anagen-to-catagen transition, signals involved in epithelial cell proliferation and differentiation are shut down (
). Since TGF-
acts in a threshold-like switch fashion, we cannot completely separate the pro-proliferation and pro-apoptosis effects of TGF-
in our model. We test four cases: 1) decreasing BMP level (BMP−) while increasing Str-TGF-
level (Str-TGF-
+), 2) BMP−/TNF+, 3) TGF-
−/TNF+, 4) TGF-
−, BMP−/TNF+ (Fig. S21). In cases 1) and 2), chances of apoptosis are almost the same as those without decreasing BMP (Figs. 4C and D, S21A and B), and it seems that the decrease of differentiation only compensates cell division and has no great effect on apoptosis. Similarly, in case 3), TGF-
−/TNF + while keeping BMP the same level, it does not change the chance of apoptosis, yet the decrease of cell division compensates the increase of differentiation (Fig. S21C). On the other hand, in case 4), with both BMP and TGF-
decreased, an increase in TNF results in a further slight increase in apoptosis to ∼0.8, in comparison with increasing TNF without decreasing BMP and TGF-
(Figs. 4D and S21D). Our simulation results imply that shutting down both the pro-proliferating and pro-differentiation signals during the apoptosis may enhance the apoptosis initiation, while shutting down either the pro-proliferating or pro-differentiation signal has little effect on apoptosis, and instead seems to compensate each other.
Testing a two-step catagen initiation hypothesis: TGF-
initiates the initial wave of apoptosis, followed by the upward propagation driven by TNF
Although experimental results demonstrate that various signals including TGF-
and TNF contribute in the initiation and propagation of the HF catagen apoptosis wave, interestingly, it is reported that TGF-
could induce catagen-like morphological changes while TNF may cause morphological changes of whole HFs that are not commonly seen in catagen HFs (
), suggesting that at least TNF is not the unique apoptosis-initiating signal in HF. Taking these experimental findings into account, and that TGF-
and TNF may have different sender/receiver cells in the HF, there is one catagen initiation hypothesis: upon the anagen-to-catagen transition, DP cells produce TGF-
—and possibly also other signals, including FGF5 (
). To test this hypothesis, we design a
cell array, where each cell carries a pBN model. BMP and TGF-
are imposed as exterior opposing signal gradients that simulate their distributions in HF Mx, and cells undergoing apoptosis produce TNF signals that are received by itself and its neighboring cells. Upon the end of anagen when BMP and TGF-
are both high, but with TGF-
below the threshold, we continue to increase TGF-
from the bottom of the domain to activate the MAPK pathway. Simulations from this pBN-equipped spatial model show that the elevation of TGF-
at the bottom of the Mx initiates an apoptosis wave (Fig. 5A), and that TNF together with strong TGF-
further drives the upward propagating apoptosis that degenerates the HF Mx (Fig. 5B and C; Video S1). Details of the spatial model can be found in Section S4. We also point out that, in reality, apoptosis causes the upward degeneration of the Mx that carries the DP, which may further contribute to the apoptosis propagation as the Str-TGF-
resource is moving up.
Figure 5The spatial simulation of apoptosis initiation and progression during catagen. During anagen, BMP and TGF-
cast opposing gradients on the Mx. Upon the anagen-to-catagen transition, TGF-
increases from the bottom of Mx, activating the Str-TGF-
input node that leads to the initiation of apoptosis (A). Apoptosis cells produce TNF, which can be received by its neighboring cells, resulting upward propagating waves of TNF and apoptosis (B and C). The strength of the signal or the probability of the apoptosis fate is shown by the color, with color bar given at the left of panel A. To see this figure in color, go online.
Predictions of gene KO and OE on HF cell fate regulations
We perform KO and OE simulations with different combinations of signal inputs, so to study the effect of each gene (node) on HF cell fate regulations, and the results are shown in Figs. 6 and S22–S43. Details of the simulation setup are provided in Section S5.
Figure 6Effects of overexpression (OE) and knockout (KO) of each gene on cell fate regulations. Each pixel shows the difference of the noted cell fate between the OE or KO simulations of the noted gene and the WT simulations, under the noted input state. Each group of simulations consists of 10,000 runs and the results are taken at the end of t = 1000. The change in the cell fate is shown by the color, with warm color (difference > 0) indicating an increase in the cell fate, and cold color (difference < 0) a decrease. (A and A′) Change in cell division due to the OE and KO of a gene. (B and B′) Change in differentiation. (C and C′) Change in apoptosis. To see this figure in color, go online.
In Fig. 6, each pixel shows the change in the noted cell fate (shown by the title) under the given signal inputs (shown by the left labels), due to the OE or KO of this gene (shown by the bottom labels), compared with the WT under the same input state. The level of cell fate change is shown by the warm/cold color indicating an increase/decrease. Fig. 6 allows us to visually identify the roles of each signal in HF cell fate regulation. In each panel, if the OE/KO mostly increases or decreases this cell fate, we highlight this gene in red/blue, respectively. If the OE/KO of a gene does not result in a clear cell fate change under most input states, we leave the gene name black. First, cell division clearly requires signals in the TGF-
-SMAD2/3-division pathway, including SMAD2/3, p-SMAD2/3, p-SMAD2/3-SMAD4 and TMEFF1. Moreover, we find that several signals in the apoptosis module also contribute to division, including c-FLIP and BCL-XL/2 in both extrinsic and intrinsic mitochondria pathway (nodes BCL-XL/2 and mBCL-XL/2). Similarly, signals in the BMP-SMAD1/5/8-differentiation pathway clearly contribute to differentiation, including SMAD1/5/8, p-SMAD1/5/8, p-SMAD1/5/8-SMAD4, and SMAD6, and c-FLIP and BCL-XL/2 from the apoptosis module also contribute to differentiation. We also notice that OE of p-SMAD2/3 and p-SMAD2/3-SMAD4 in the TGF-
-SMAD2/3-division pathway may decrease the chances of both differentiation and apoptosis, while the KO of these signals appears to increase the chance of differentiation more than that of apoptosis. A similar effect is observed in p-SMAD1/5/8 and p-SMAD1/5/8-SMAD4 from the BMP-SMAD1/5/8-differentiation pathway, that OE of these signals decreases the chances of both division and apoptosis, while the KO of these signals increases the chance of division more than apoptosis.
As for apoptosis, we identify several signals/nodes in the apoptosis module that contribute to apoptosis, including TRADD, TRAF/RIP, FADD, CASP8/10, CASP3/7, BID, JNK, DAXX, p53, BAX/BAK, CASP9, and MOMP. As a connector between the division/differentiation and apoptosis modules, SMAD7 also contributes to apoptosis while inhibiting division and differentiation. Effects of several signals in the TGF-
-SMAD2/3-division and BMP-SMAD1/5/8-differentiation pathways have been discussed above. Other signals in the apoptosis module can be divided into three groups based on their effects on cell fate regulation. The first group includes c-FLIP and BCL-XL/2, which promote both division and differentiation, and inhibit apoptosis, as discussed above. The second group includes MEK1/2, RAS, and RAF. KO of these signals slightly increases the chance of apoptosis, although OE of the signals shows no clear effect in changing the chance of apoptosis. This demonstrates that these signals act slightly in an apoptosis-inhibiting fashion. Note that these signals also inhibit both division and differentiation, and we find that, surprisingly, they act in a way that inhibits all cell fates; in other words, KO of these signals may increase all three cell fates. Considering the modeled undetermined fate in our model, it might be possible that some circuits involving these signals are currently absent, which should be able to provide the cell the missing information in cell fate decision. We suggest future experiments to investigate in this direction. Finally, there is a third group in the apoptosis module, including IKK and NF-κB. While OE of IKK and NF-κB dominantly inhibit apoptosis as well as division and differentiation, KO of them slightly increases apoptosis under some input states, while slightly decreasing apoptosis under some other input states. We consider that IKK and NF-κB have a hybrid role in regulating apoptosis, and we highlight them in purple in the KO plot of apoptosis (Fig. 6C′).
In the following, we briefly compare our simulation findings with known experimental results from the literature.
SMAD6, TMEFF1
The BMP-SMAD1/5/8-differentiation and TGF-
-SMAD2/3-division pathways inhibit each other through SMAD6 and TMEFF1. The inhibition of TMEFF1 from TGF-
-SMAD2/3 to BMP-SMAD1/5/8 pathway in HF stem cell cells is reported in (
). On the other hand, the regulating role of SMAD6 in HF is unclear, and we would point out that there is evidence of SMAD6 inhibiting on the BMP-SMAD1/5/8 pathway in other systems (
). However, our simulations predict that, without such an inhibitory node from BMP to the TGF-
pathway, cells will dominantly commit cell cycling instead of differentiation. This indicates that even the regulatory role of SMAD6 still needs to be confirmed, and there should be a BMP to TGF-
inhibitory node that serves the regulation in the same way.
The two apoptosis pathways
Str-TGF-
may regulate cell apoptosis through two pathways: Str-TGF-
-RAS-MEK1/2, or Str-TGF-
-DAXX-JNK. As discussed above, KO any of RAS, RAF, or MEK1/2 slightly increases all chances of division, differentiation, and apoptosis. Such a dual role results mainly from the dual role of NF-κB, which can be anti-apoptotic through activating BCL-XL/2 as well as the activation of self-inhibition possibly through c-FLIP (
cFLIPL inhibits tumor necrosis factor-related apoptosis-inducing ligand-mediated NF-κB activation at the death-inducing signaling complex in human keratinocytes.
), and meanwhile it also activates SMAD7, which inhibits both the BMP-SMAD1/5/8-differentiation and TGF-
-SMAD2/3-division pathways. Next, for the Str-TGF-
-DAXX-JNK pathway, we point out that, although the (Str-)TGF-
-JNK pathway is more established in literature, the role of DAXX is still under debate, and there is evidence on both the activating and inhibiting roles of DAXX on cellular apoptosis, through the interaction with the JNK pathway (
SMA7 is activated by NF-κB and can shut down both the BMP-SMAD1/5/8 and TGF-
-SMAD2/3 pathways. In addition, it may also have an inhibitory effect on IKK, which further inhibits NF-κB (Fig. 1B). Simulation results show that the inhibitory regulations of SMAD7 on the BMP-SMAD1/5/8 and TGF-
-SMAD2/3 pathways are more prominent in comparison with its inhibitory role in the apoptosis pathway. It is reported that OE of SMAD7 may result in delayed and aberrant HF formation, implying its inhibitory effect on cell division (
), supporting its dominant pro-apoptosis effect shown by our simulations. BCL-XL/2 are known as apoptosis inhibitors, and our simulations validate their apoptosis-inhibiting effect. Experimental evidence demonstrates that BCL-2
null mice showed a retarted anagen development, and hair shafts were produced 1.2–1.4 times more slowly compared with WT mice (
bcl-2 deficiency in mice leads to pleiotropic abnormalities: accelerated lymphoid cell death in thymus and spleen, polycystic kidney, hair hypopigmentation, and distorted small intestine.
). However, we point that other experiments on BCL-2 OE show a seemingly dual role on apoptosis regulation: while BCL-2 OE is reported to reduce ultraviolet B-induced apoptosis in mutant mice, BCL-XL/2 OE may accelerate catagen development, possibly through the interaction with the FGF5 pathway (
In recent years, HFs have become a popular model system in the study of morphogenesis, growth, and stem cell biology. While there have been many modeling works on the morphogenesis and growth control of HFs, modeling studies on the HF gene regulation mechanisms and the cell fate control have been limited. In this work, we developed a Boolean model from the literature, then used scRNA-seq data to refine it to a probabilistic Boolean model for HF epithelial cells. We used the pBN to investigate the cell fate regulatory mechanisms in HF epithelial cells, focusing on the dual roles of TGF-
and its interaction with other signals. We propose a threshold-based switch on TGF-
regulation on HF cell fate, and simulation results show qualitative match with experimental results (
). Further attractor and stable motif analysis on the associated Boolean model show that the division and differentiation cell fates result from fixed-point steady states, while apoptosis leads to complex attractors. Then, we use the model to investigate the signaling mechanisms behind two critical “checkpoints” in the HF growth cycle: the telogen-to-anagen transition and the anagen-to-catagen transition. At the telogen-to-anagen transition, it is known that a typical high inhibitor/low activator profile during the refractory telogen cannot enable the telogen-to-anagen transition. Instead, the HF needs to possess the low inhibitor/low activator profile during competent telogen. Our simulations show concordance with this telogen-to-anagen transition mechanism, and reveal its dependence on the cross-inhibition between the BMP-SMAD1/5/8-differentiation and TGF-
-SMAD2/3-division pathways through the two critical regulatory nodes TMEFF1 and SMAD6. At the anagen-to-catagen transition, our simulations on a single cell show that both TGF-
and TNF can initiate apoptosis, with TNF being more efficient. Two other hypotheses related to the anagen-to-catagen transition are also investigated: first, shut-down of both division and differentiation signals may further enhance apoptosis; next, the apoptosis wave may be first initiated by TGF-
, then propagate through TNF produced by cells undergoing apoptosis as well as TGF-
. Finally, we perform KO/OE simulations of the genes, and compare the results with some experimental results from the literature. In particular, we find that signals in the apoptosis module may enhance apoptosis (TRADD, TRAF/RIP, FADD, CASP8/10, CASP3/7, BID, JNK, DAXX, p53, BAX/BAK, CASP9 and MOMP, and SMAD7), or inhibit apoptosis (c-FLIP and BCL-XL/2), or have hybrid effects on cell fate regulation (MEK1/2, RAS, RAF, IKK, and NF-κB).
While our pBN model offers a useful approach in the study of HF biology and provides interesting insights to the HF epithelial cell fate regulation dynamics, there are a few directions that can be further improved in the future. Currently the GRN in our model only involves three signals, TGF-
, BMP, and TNF, and together they form a minimal set that allows cell to undergo division, differentiation, and apoptosis. In fact, cell fate regulatory mechanisms in the HF are highly complicated and involve many other signals, and together they drive the cell fate decisions at different stages of an HF growth cycle. To name a few, it is well-known that WNT is important in initiating and maintaining the anagen, sonic hedgehog (SHH) induces HF formation during morphogenesis and it also guarantees the formation of HF ORS during anagens, and FGF5 is a signature signal in catagen initiation. An important and urgent next step is to extend the current three-signal GRN to include at least these well-established HF cell fate regulatory signals. This may help us better understand the nature of the HF GRN nonlinear dynamics. Moreover, as a mini-organ, an HF consists of different compartments, and cells of different types actively “talk” to each other through different signals, so to orchestrate the HF growth and maintain its functions. Recently, new methods have been developed to infer the cell’s communication events (
). A combined model of both inter-cellular signaling communications and intra-cellular gene regulations might help people better understand how the HF cell fate decisions are coordinated at both the cellular and the whole-organ level.
) from HF Mx epithelial cells to refine the Boolean network. While apoptosis is first induced in Mx cells upon the anagen-to-catagen transition, during late telogen, TGF-
acts on HF stem cells to initiate anagen. Therefore, to better understand the signaling events behind HF cell fate decisions and the HF cyclic growth dynamics, we need to investigate the HF stem cells and their direct progeny ORS cells, as well as data from telogen HFs (
). On the model refinement side, in recent years, various GRN inference methods based on single-cell data have been developed, including several ones that focus particularly on the inference of Boolean network models (
). In a future extension of this work, using multiple inference methods and comparing their results may improve the accuracy of the refined model. We also point out that many of these Boolean model inference methods were originally developed for qPCR data rather than scRNA-seq data, which presents a challenge to the inference accuracy due to the high level of noise typically presented in the scRNA-seq data. Applying imputation methods to the scRNA-seq data before applying the inference method might partially rescue it. In addition, we note that the PB scores and the probabilities in the pBN do not reflect true biological meaning (see discussion in Section S3). While our data-refined pBN modeling approach provides insights to the underlying mechanisms of cell fate decisions and useful predictions, further functional experiments are needed to validate the GRN. Finally, the Boolean model assumes each gene (node) to have two discrete states, ON or OFF, which oversimplifies the gene regulation dynamics. Thus, in the future, we should develop continuous models based on the GRN structure developed in this paper. The ordinary differential equation (ODE) type of GRN models might be a good option in this direction, with recently available ODE-based single-cell data-inference methods (
Q.W. conceived the original idea and designed the model. K.D. and Q.W. analyzed the data, performed the numerical simulations, and wrote the article.
Acknowledgments
The authors thank Dr. Simon Joost and Dr. Maria Kasper for providing their original data and their results of diffusion map and pseudotime analysis, Dr. Nicola D. Wilson and Dr. Berthold Göttgens for their help with the pseudotime Boolean method, and Dr. Maksim Plikus for valuable discussions on the modeling development and results. The authors acknowledge support from the National Science Foundation grant DMS-1951184 to Q.W. Computations were performed using the computer clusters and data storage resources of the HPCC at UC Riverside, which were funded by grants from NSF (MRI- 1429826) and NIH (1S10OD016290-01A1).
Network modeling of TGFβ signaling in hepatocellular carcinoma epithelial-to-mesenchymal transition reveals joint sonic hedgehog and Wnt pathway activation.
Boolean model of growth signaling, cell cycle and apoptosis predicts the molecular mechanism of aberrant cell cycle progression driven by hyperactive PI3K.
Smad7 sensitizes tumor necrosis factor–induced apoptosis through the inhibition of antiapoptotic gene expression by suppressing activation of the nuclear factor-κB pathway.
cFLIPL inhibits tumor necrosis factor-related apoptosis-inducing ligand-mediated NF-κB activation at the death-inducing signaling complex in human keratinocytes.
bcl-2 deficiency in mice leads to pleiotropic abnormalities: accelerated lymphoid cell death in thymus and spleen, polycystic kidney, hair hypopigmentation, and distorted small intestine.
00428-3&id=gr1.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00428-3&id=gr2.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00428-3&id=gr3.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00428-3&id=gr4.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00428-3&id=gr5.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}
00428-3&id=gr6.jpg&_x_tr_sl=en&_x_tr_tl=fr&_x_tr_hl=fr&_x_tr_pto=wapp){kind=link}