Hepatocytes have genetic markers that influence EMT
In order to analyzing dynamical properties of EMT in HCC context, we first searched for genetic markers and molecular features that distinguish hepatocytes from other epithelial cells. We then incorporated such interactions to a previously published and validated model of EMT proposed by Méndez-López et al37. In this model we added specific transcription factors (TFs) that have been described in hepatocyte differentiation such as HNF4A, HNF1A, FOXA2, and HNF6. These genetic markers belong to a cross-regulatory network that controls the development and adult function of the liver38,39. Particularly, it has been reported that HNF4 and HNF1A are negative regulators of Snail and Slug40. Moreover, we incorporated TFs and interactions that are related to stemness and endodermal differentiation. In this sense, we added to the network OCT4, SOX2, NANOG and KLF4 because they are central regulators for the induction and maintenance of stem cells41–43. It is noteworthy to mention that such TFs are involved in the emergence of a small stem-like cell population present in injured liver of mice44.
Other genetic markers that have been incorporated to the network are GATA6, SOX9, b-catenin, YAP1 as well as miR-200a,b,c, and miR34a. Concerning each marker, GATA6 is essential for visceral endoderm differentiation and embryonic development of the liver 45,46. SOX9 expression has been associated with proliferation and stem cell features in HCC47. SOX9 is expressed in liver stem/progenitor cells but not in hepatocytes48,49. Wnt/b-catenin signaling pathway plays a pivotal role in liver development and regeneration50. In particular, b-catenin regulates EMT in HCC cells through double-negative feedback with HNF4A51. Regarding YAP1, it forms a circuit with SNAI1 and HNF4A52, to regulate EMT53, proliferation, and differentiation in HCC cells54. The miR-200 family (miR-200a,b,c) and miR34-a are involved in EMT55, reprogramming56, and differentiation57. MiR-34a plays a role in EMT58. SNAI1 and HNF4A regulate the expression of the miR-200 family and miR-34a. Grounded on experimental evidence of a total number of 240 papers (Additional File 1), we were able to postulate a regulatory network of 43 nodes that represents hEMT in the context of HCC (Fig. 1). We use all the information presented in this network to build a Boolean model of 43 logic rules, available in additional file 2 (Additional File 2).
hEMT regulatory network is a source of phenotypic diversity
The network presented above, includes much of the complexity associated with hEMT. However, we seek to understand what is the effect of mutations observed in cancer cells on hEMT. For this reason, we decided to find the necessary and sufficient nodes to control the global dynamics of the GRN. To do this, we proceed to compact the network model using the algorithm proposed by Véliz-Cuba 59. Such procedure uses discrete math operations to simplify linear interactions, conserving all non-linear motives that contribute to the system dynamics, like positive and negative feedback loops (see Methods). Using this algorithm, we obtained a reduced network of 23 nodes (Fig. 2). In the same way, we validated our reduction by using the GINsim software 60. All Boolean logic rules of this reduced network are presented in additional file 3 (Additional File 3).
After reduction, we needed to know whether the simplified model of the hEMT network had preserved the essential dynamic information of the full model of GRN. To this end, we identified the attractors of both Boolean models using the BoolNet R package. In particular, for the simulation of the GRN full model, 5,000,000 configurations of the network as initial conditions were randomly selected, and it was found that the model converges to only eight attractors that correspond to the following phenotypes: Senescent hepatocytes (SH), senescent mesenchymal cells (SM), quiescent hepatocytes (QH), quiescent mesenchymal cells (QM), proliferative hepatocytes (PH), proliferative mesenchymal cells (PM), proliferative stem-like hepatocytes (PSH), and proliferative stem-like mesenchymal cells (PMS) (Fig. 3a).
The mesenchymal phenotype was identified by the activation of SNAI1, SNAI2, TWIST1, ZEB1, FOXC2 15,61. The activation of HNF4A, HNF1A, HNF6 and FOXA2 correspond to the hepatocyte phenotype 38,62. In the same way, the senescent phenotype was identified by the activation of p53, p16, p21 and RB 63,64; proliferation was identified by E2F1 activation, cyclin D, and RB inactivation 65,66; quiescence was identified by RB and p21 activation, as well as by the absence of proliferation markers like E2F1, cyclin D, and the absence of senescence markers such as p53 and p16 67. Finally, the stem-like phenotype was associated with the activation of OCT4, SOX2, and NANOG 43,68.
As expected, the reduced model presented equivalent attractors (Fig. 3b), which indicates that set of nodes of this model is sufficient to describe the network dynamics. Subsequently, we compared the size of the basin of attraction for each phenotype. This was done for both the full model (Fig. 3c) and the reduced model (Fig. 3d). We observed that such results are qualitatively congruent with each other. Remarkably, PMS and PSH presents stem-like features given by the activity of NANOG, OCT4 and SOX2 (Fig. 3a and 3b), which suggests that such phenotypes have self-renewal properties. Moreover, the size of the basin of attraction of PMS and PSH indicates that such phenotypes may be less frequent than others (Fig. 3c and 3d), which is congruent with the low prevalence of stem cells in nature 69. Collectively, these results indicate that the reduced model captures the biological richness of the extended model. Therefore, it could be used to assess the effect of mutations on hEMT plasticity.
Phenotypes produced by hEMT GRN are robust and physiologically feasible
The results presented above show that hEMT GRN reproduces a set of well-defined hepatic phenotypes. In a physiological context, phenotypes persist despite slight molecular changes produced by viral infections, thermal changes, and osmotic pressure variations, among other stimuli 70. In other words, phenotypes are robust. For this reason, we determined whether these phenotypes produced by the hEMT model could persist in the presence of perturbations. To do this, we tested the robustness of the GRN in presence of fluctuations by starting Boolean simulations with all states of the network as an initial conditions, and then applying stochastic noise to the logical rules of the nodes in order to determine whether the state belongs to the same basin of attraction or whether it transits to a different basin 71. As a result of this procedure, we determine the probabilities that an hEMT-produced phenotype remains fixed or to change towards other phenotypes randomly (Fig. 4a). Numerically, we found the probabilities that the phenotypes are maintained in the face of fluctuations are higher than 0.9, which indicates that such phenotypes are robust and, therefore, physiologically feasible (Fig. 4b).
Experimental observations on hEMT are mechanistically explained by the GRN
After determining the feasibility of the minimal hEMT model under stochastic fluctuations, we test its effectiveness as a predictive tool. To this end, we simulate a series of loss-or-gain-of-function mutations on a particular set of genes like SNAI1, SNAI2, HNF1A and HNF4A as well as oncogenes as YAP1 and b-catenin. In the same way, we studied the effects of targeting tumor suppressors such as p53, p21, p16 and RB (Fig. 5a and Additional File 4). To simulate the gain of function, we set the value of the target node as one. Similarly, to simulate the loss of function, we set the value of the target node to zero for all time-steps of simulation. Next, we identified the attractors and basin sizes of each mutated network. Finally, we analyze the impact (\(\varDelta {S}_{i}\)) that each mutation has on the WT attractor landscape by subtracting the value of each wild-type basin (\({S}_{WT}\)) from the value of the mutated basin (\({S}_{i}\)), that is: \(\varDelta {S}_{i}={S}_{i}-{S}_{WT}\) (Fig. 5b).
As a result of these assays, the model was able to qualitatively reproduce several experimental observations, such as either the inhibition of HNF4A and HNF1A or the overexpression of SNAI1 and SNAI2 increase the mesenchymal phenotype (Fig. 5b and Table 1). The opposite occurs when HNF4A and HNF1A are overexpressed, increasing the hepatocyte phenotype (Fig. 4d and Table 1). In the same direction, inhibition of p53, p16 and p21 reduces senescent phenotypes (Fig. 5b and Table 1), while the opposite occurs when p53 and p16 are overexpressed (Fig. 5b and Table 1). Furthermore, the model showed that proliferation can be triggered when p21, p16, p53, and RB are inhibited or when YAP1 and b-catenin are overexpressed (Fig. 5b and Table 1). Finally, the model also showed that the constitutive activation of YAP1 increases the stem-like phenotype (Fig. 5b and Table 1). Collectively, these results validate our qualitative model as a predictive tool.
Table 1
Validation of the hEMT model
Mutation* | Prediction** | Experimental outcomes | References |
KO HNF4A | This condition increases mesenchymal cells. | HNF4A silencing increased migratory capacity and the expression of mesenchymal markers. | 72 |
KO HNF1A | This condition increases mesenchymal cells. | This increases cell’s migratory capacity and the expression of mesenchymal markers. | 72 |
KO SNAI1 | This condition reduces mesenchymal cells. | SNAI1 silencing decreased migratory capacity and the expression of mesenchymal markers. | 27,73 |
KO p53 | This condition reduces senescence. | p53 deletion reduces senescence. | 63,74,75 |
KO p16 | This condition reduces senescence. | p16 deletion reduces senescence. | 76 |
KO p21 | It reduces senescence and increases proliferation. | p21 deletion reduces senescence and increases proliferation. | 77,78 |
KO RB | This condition increases proliferation. | RB deletion increases proliferation. | 77,79,80 |
OE HNF4A | This condition reduces mesenchymal cells. | HNF4A overexpression increases epithelial morphology and reduces motility as well as invasive capacity. | 27,72 |
OE SNAI1 | This condition increases mesenchymal cells. | SNAI1 overexpression increase mesenchymal phenotype. | 27,73,81 |
OE SNAI2 | This condition increases mesenchymal cells. | SNAI2 overexpression increase mesenchymal phenotype. | 40,82 |
OE b-catenin | This condition increases proliferation. | b-catenin activation promotes proliferation. | 83,84. |
OE YAP1 | This condition increases stem-like phenotype. | YAP1 activation increase stem-like phenotype. | 85 |
OE p53 | This condition increases senescence. | p53 activation promotes senescence. | 86 |
OE p16 | This condition increases senescence. | p16 overexpression promotes senescence. | 76 |
*KO means ‘knockout’ and OE is ‘overexpression’. **All predictions were taken from Fig. 5b. |
Hepatocytes and other epithelial cells may be similarly transformed into CSCs
The analysis of mutations on hEMT GRN showed that there are genes capable of biasing the process towards proliferative phenotypes, including the PMS phenotype (Fig. 6a and 6b). Consequently, we decided to investigate the effect on the robustness of attractors produced by the loss of function of the tumor suppressors p53, p21, p21 and RB, as well as the overexpression of YAP1 and b-catenin oncogenes. Regarding the phenotype robustness, we observed that all simulated mutations enhance proliferative phenotypes, including PMS (Fig. 6c). Conversely, senescent phenotypes were predominantly affected by p53 inactivation (Fig. 6d). The decreased robustness of the senescent phenotypes by the p53 mutation may explain the senescence escape experimentally observed 11. Interestingly, mutations on tumor suppressors as well as oncogenes favor the appearance of proliferative stem-like phenotypes (Fig. 6d). We hypothesize that this augment in stem-like phenotypes corresponds to the appearance of CSCs, represented by PMS and PSH phenotypes.
It is not clear whether this property is exclusive to hepatocytes or not. To explore this issue, we compared the impact of each mutation (\(\varDelta {S}_{i}\)) on the WT attractor landscape of hEMT model against Méndez-López et al 37 observations in generic EMT model (Fig. 6e). Qualitatively, the behavior of hepatocytes was similar to other epithelial cells. However, it is notable that hepatocytes were less sensitive to SNAI2 mutations. In fact, SNAI2 did not induce severe alterations in senescent phenotypes from the hEMT model compared to other epithelial cells (Fig. 6e). Furthermore, we noted that KO p53 strongly reduced senescent phenotypes in liver compared to other epithelial cells (Fig. 6e). Collectively, these results suggest that there might be a common mechanism in the appearance of CSCs in epithelial cells. However, CSC emergence may be affected by the phenotype-specific regulators of each cell type.
Mutations in oncogenes and tumor suppressor genes affect EMT plasticity and generate CSCs
Next, we determine what the most likely phenotypic transitions are. To investigate this point, we use the mean first passage time (MFPT), which is a metric to determine how difficult it is for a GRN to take the first step in a differentiation trajectory towards another phenotype. In general, the MFPT is used to calculate the net transition rate, which determines how easy it is to move from one state to another (see Methods). As a result of this procedure, we observed no transitions from senescent phenotypes to proliferative or quiescent phenotypes in WT model (Fig. 7a), which confirms the stability of the senescent phenotypes suggested by the size of the attraction basins (Fig. 3c and 3d). We also noted that transitions towards proliferative and stem-like phenotypes are relatively low (Fig. 7a). In the WT model, the attractors associated with the quiescent and senescent phenotypes are the most stable.
Finally, we examined how mutations could affect the probabilities of hEMT. To do this, we simulated common alterations in HCC and other cancers like loss-of-function mutations of RB, p21, p16, and p53 and constitutive activation of b-catenin, YAP1, NF-kB, and TGF-b. As a result of these simulations, we found that aberrant activation of b-catenin, YAP1, NF-kB, and TGF-b increased the index of net transition of proliferative hepatocytes to proliferative mesenchymal cells (Fig. 7b), in concordance with previous experimental observations 16,28,31. On the other hand, the loss-of-function of tumor suppressor genes like p16, p21 and p53 also increased this transition (Fig. 7b), in agreement with in vitro assays 22,24,87. All simulated alterations promoted the transition from hepatocytes with stem-like features to mesenchymal stem-like phenotype. However, only p53 loss-of-function mutation and constitutive activation of NF-kB and TGF-b favored the transition of proliferative hepatocytes to proliferative mesenchymal stem-like phenotype (Fig. 7b). Thus, these mutations may be essential to increase the prevalence of CSCs.