{'fraction_non_alphanumeric': 0.033516189456783516, 'fraction_numerical': 0.018196414236018196, 'mean_word_length': 5.296686324066274, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics of hepatocyte-cholangiocyte cell-fate decisions during liver development and regeneration

Summary The immense regenerative potential of the liver is attributed to the ability of its two key cell types – hepatocytes and cholangiocytes – to trans-differentiate to one another either directly or through intermediate progenitor states. However, the dynamic features of decision-making between these cell-fates during liver development and regeneration remains elusive. Here, we identify a core gene regulatory network comprising c/EBPα, TGFBR2, and SOX9 which is multistable in nature, enabling three distinct cell states – hepatocytes, cholangiocytes, and liver progenitor cells (hepatoblasts/oval cells) – and stochastic switching among them. Predicted expression signature for these three states are validated through multiple bulk and single-cell transcriptomic datasets collected across developmental stages and injury-induced liver repair. This network can also explain the experimentally observed spatial organization of phenotypes in liver parenchyma and predict strategies for efficient cellular reprogramming. Our analysis elucidates how the emergent dynamics of underlying regulatory networks drive diverse cell-fate decisions in liver development and regeneration.


INTRODUCTION
The liver, the largest internal organ in the body, performs key physiological functions. It possesses remarkable regenerative ability and is capable of restoring its mass, architecture, and function completely after injury (Gadd et al., 2020;Kopp et al., 2016). Both the major cell types seen in liver parenchyma -hepatocytes and cholangiocytes (biliary epithelial cells) -have been shown to be capable of dividing extensively (Kopp et al., 2016) and trans-differentiate into one another, thus contributing to liver regeneration (Deng et al., 2018;Schaub et al., 2018;Yanger et al., 2013).
In the liver, hepatocytes and cholangiocytes perform very different functions. Although the former executes most metabolic functions, including bile secretion; the latter are biliary epithelial cells that line the bile duct tubules and serve to transport bile from liver to the small intestine. Developmentally, hepatocytes and cholangiocytes are both formed from a common progenitor cell type known as hepatoblasts. Hepatoblasts can co-express markers for both hepatocytes (HNF4a, CK18) and cholangiocytes (CK19) (Gordillo et al., 2015), a characteristic trait of such bipotent cells witnessed during embryonic development (Zhou and Huang, 2011). Cellular plasticity seen in vivo between cholangiocytes and hepatocytes suggests that although they may be terminally differentiated, they retain the ability to trans-differentiate into one another during injury-induced repair (Deng et al., 2018;Schaub et al., 2018). These observations suggest that these differentiated cells in liver can carry permissive chromatin from their progenitors, which may be important for their developmental and reprogramming competence (Li et al., 2020). However, the intracellular and tissue-level dynamics of cell-state transitions between these cell types in a liver remains unclear.
Phenotypic plasticity is a hallmark of injury response too. For instance, during chronic liver injury, a subset of resident liver cells are capable of undergoing epithelial-mesenchymal transition (EMT) and later regenerate hepatocytes through mesenchymal-epithelial transition (MET) (Choi and Diehl, 2009). EMT/MET has been well-investigated in such pathophysiological conditions (Nieto et al., 2016). However, in another manifestation of plasticity, the cholangiocytes can transition to hepatocytes via an intermediate biphenotypic state which expresses both hepatocytic and biliary markers (Deng et al., 2018;Manco et al., 2019), as identified via lineage tracing. Importantly, these reprogrammed cells are more proliferative and less apoptotic than the native hepatocytes (Manco et al., 2019). The biphenotypic state seen is reminiscent of earlier observations about the existence of oval cells -bipotent cells that expressed hepatoblast marker AFP and could

RESULTS
Multistability in gene regulatory network underlying hepatocyte-cholangiocyte cell-fate commitment and plasticity First, we identified a core gene regulatory network that is integral to the observed plasticity among hepatocytes, hepatoblasts and cholangiocytes. Various molecular players have been implicated in remarkable cellular plasticity of liver cells both in the context of liver development (Gé rard et al., 2017;Lau et al., 2018;Zong and Stanger, 2011) and regeneration (Li et al., 2020). Among these, TGFb signaling is of key importance in the differentiation of cholangiocytes from the hepatoblasts. High levels of TGFb signaling is required near the portal vein for differentiation of biliary cells in vivo (Clotman et al., 2005). Consistently, TGFb can mediate hepatocyte to cholangiocyte trans-differentiation in vivo (Schaub et al., 2018), and can suppress the transcription and activity of HNF4a, a well-known hepatocyte inducer (Cozzolino et al., 2013;Li et al., 2000;Lucas et al., 2004). Overexpression of TGFBR2 in hepatoblasts in vitro can drive them into a cholangiocyte fate (Takayama et al., 2014). Thus, TGFb signaling can be thought of as a 'master regulator' of cholangiocyte cell-fate. Similarly, a 'master regulator' of hepatocyte cell fate is c/EBPa, whose overexpression transcriptionally inhibits TGFBR2 and can drive hepatocyte differentiation in hepatoblasts (Takayama et al., 2014). Furthermore, suppression of c/EBPa can stimulate biliary cell differentiation via increased Hnf6 and Hnf1b expression in periportal hepatoblasts (Yamasaki et al., 2006). In these hepatoblasts, overexpression of TGFBR2 inhibits c/EBPa, thus forming a 'toggle switch' or mutually inhibitory feedback loop between TGFBR2 and c/EBPa. Such 'toggle switches' are hallmarks of cellular decision-making between various 'sibling' cell fates (Zhou and Huang, 2011).
Another important player in biliary development is SOX9, which can repress the expression of both c/EBPa and TGFBR2 in mature biliary cells (Antoniou et al., 2009;O'Neill et al., 2014). SOX9, together with SOX4, can coordinate development of the bile duct (Poncy et al., 2015), during which it is upregulated by TGFb signaling through the Jagged1-Notch axis (Wang et al., 2018). SOX9, similar to c/EBPa and TGFb/TGFBR2 signaling, has been proposed to self-activate either directly or indirectly (Du et al., 2018;Duan and Derynck, 2019;Friedman, 2015). Put together, these interactions constitute a gene regulatory network comprising SOX9, c/EBPa and TGFBR2 ( Figure 1A) which is involved in hepatocyte-cholangiocyte cell-fate decisions in the liver. Our goal here is not to identify a comprehensive network, but a minimal network motif that can potentially explain diverse instances of cell differentiation and reprogramming seen in liver ll OPEN ACCESS 2 iScience 25, 104955, September 16, 2022 iScience Article development and injury repair. Thus, the nodes represented here -SOX9, c/EBPa and TGFBR2 -can be viewed as proxies for their corresponding co-factors and their regulons implicated in cell-fate commitment.
To understand the distinct phenotypes enabled by the gene regulatory network involving c/EBPa, TGFBR2 and SOX9 ( Figure 1A), we performed dynamical simulations on this regulatory network, using an ensemble of kinetic parameter sets for a set of ordinary differential equations (ODEs) representing the regulatory interactions within a network. The parameter sets sampled randomly from a biologically plausible regime capture the inherent variability in a cell population. For each parameter set, various possible steady-state (phenotype) combinations are collated to identify the robust dynamical features emerging from this network topology (see STAR methods).
For the ensemble of steady-state solutions thus obtained, hierarchical clustering suggested the evidence of three distinct cell-states, through a heatmap Figure 1B). A c/EBPa-high state can be attributed to a  iScience Article mature hepatocyte (Akai et al., 2014), whereas a SOX9-high state can be mapped on to a cholangiocyte profile (Dianat et al., 2014). Thus, progenitor-like (hepatoblasts or oval cells) state is likely to have low levels of both c/EBPa and SOX9, with possible enrichment of TGFBR2 (Takayama et al., 2014).Further, we observed a continuum between the progenitor-like and cholangiocyte phenotypes, pointing toward a gradient of cholangiocyte fate determination and maturation with possible intermediate states, as experimentally reported (Yang et al., 2017). Further credence for three states was obtained via principal component analysis (PCA), through which three distinct clusters of points were visually discernible ( Figure 1C). The c/EBPa-high state (cluster I) had relatively low levels of SOX9 and TGFBR2. The TGFBR2-high state (cluster II) had moderate to low levels of SOX9 and low levels of c/EBPa. The SOX9-high state (cluster III) had low levels of both TGFBR2 and c/EBPa ( Figure 1D). These observations were corroborated by experimental data from transcriptomic analysis of hepatocytes and cholangiocytes (GSE28891) (Shin et al., 2011) ( Figure 1E). As predicted by our model, levels of c/EBPa were significantly higher in hepatocytes than in cholangiocytes, but those of SOX9, TGFb and TGFBR2 were significantly higher in cholangiocytes than in hepatocytes ( Figure 1E). In another dataset where the adult hepatocyte and cholangiocyte signatures were found to be significantly enriched in hepatocytes and cholangiocytes respectively, the c/EBPa regulon activity was higher in hepatocytes than in the biliary cell population whereas the SOX9 regulon and overall TGFb signaling followed the opposite trend (GSE20498; Figure 1F) (Cullen et al., 2010). Finally, we show that in yet another dataset consisting of hepatocytes and cholangiocytes (Raven et al., 2017), c/EBPa levels were significantly higher in hepatocytes than in cholangiocytes whereas TGBFR2 expression, TGFB1 expression, and SOX9 regulon activity were higher in biliary cells compared to hepatocytes (GSE98034; Figure S1A). This analysis suggests that multistable dynamics of this gene regulatory network investigated here can allow for three distinct phenotypes -hepatocytes, cholangiocytes and bipotent progenitors (hepatoblasts/oval cells) -whose molecular footprints as predicted here are corroborated by diverse high-throughput transcriptomic datasets.
Next, we quantified pairwise correlation between these master regulators, and noted that while c/EBPa is negatively correlated with SOX9 and TGFBR2, SOX9 and TGFBR2 correlate positively ( Figure S1B) in our simulations. This trend is consistent with observations in a population of primary hepatocytes and hepatocyte progenitor cells ( Figure 1G; GSE105019) (Fu et al., 2019) and with the evidence during the hepatocyteductal trans-differentiation (O'Neill et al., 2014). Further, we plotted histograms for gene expression values of c/EBPa, SOX9 and TGFBR2 obtained from simulations, revealing their bimodality ( Figure 1H). Such bimodality was validated from single cell RNA-seq performed specifically on a population consisting of hepatocytes, hepatoblasts and cholangiocytes ( Figure 1I; GSE90047) (Yang et al., 2017).
The above-mentioned observations are largely also seen in simulations of an expanded gene regulatory network that includes a separate node showing the TGFb ligand itself (Figures S1C-S1E). Besides allowing for the three above-mentioned phenotypes, this extended network also enabled another one with high levels of c/EBPa and SOX9 ( Figure S1E). This phenotype may correspond to bipotent SOX9 + hepatocytes as has been observed in various context of liver injury and regeneration (Han et al., 2019;Tanimizu et al., 2014). Put together, we demonstrate that the network motif involving c/EBPa, SOX9 and TGFBR2 can recapitulate various cell phenotypes observed in liver development and regeneration: hepatocytes, cholangiocytes and their bipotent progenitors.
Consistent role of c/EBPa, SOX9, TGFb signaling during cell fate decisions in liver development, regeneration and reprogramming To gain further confidence in the network curated and our model predictions, we quantified the activity of various metrics relevant to our curated gene regulatory network in a single-cell RNA-seq dataset (GSE90047) for sorted hepatoblasts, hepatocytes and cholangiocytes from E10.5-E17.5 mouse fetal livers (Yang et al., 2017). As cells differentiated from hepatoblasts to hepatocytes, the adult hepatocyte signature and the c/EBPa regulon activity increased, whereas TGFb signaling levels consistently decreased (Figure 2A). Similarly, as cells differentiated into cholangiocytes, activity for adult biliary signature and SOX9 regulon increase ( Figure 2B). TGFb signaling did not show an expected decrease; instead its variance increased. In the extended network, nonetheless, mature cholangiocytes had visibly lower levels of TGFb than the bulk of developing cholangiocytes (Figures S1C and S1E), as observed in matured cholangiocytes (Antoniou et al., 2009). These observations were further corroborated by plotted the cell-fate trajectories at an individual cell level on a 3D PCA space, clearly showing the increase of c/EBPa levels in the hepatocyte branch and that of SOX9 higher in cholangiocyte branch respectively ( Figure 2C). However,  Figure S2C). Conversely, a subset of adult stem cells had high TGFb signaling and SOX9 regulon activities, indicating a possible commitment toward the cholangiocyte fate ( Figure S2C).
It was been earlier reported that progenitor cells generally have a larger transcriptional diversity quantified at a single-cell level (Gulati et al., 2020). To assess the transcriptional diversity in the cells, we computed the Shannon entropy values for individual cells (Teschendorff and Enver, 2017). Intriguingly, we found out that entropy is the highest in the context of the progenitor cell and significantly lower in the hepatocytes indicative of a more differentiated state ( Figure 2D).
After investigating developmental datasets, we investigated samples from alcoholic steatosis liver, where we observed that the steatosis liver samples showed a strong enrichment of the adult hepatocyte signature along with high c/EBPa expression (GSE102683; Figure 3A). This observation is expected as hepatocytes form a majority of the liver mass. EpCAM has been reported to be expressed in both mouse normal cholangiocytes and oval cells. TROP2, is expressed exclusively in oval cells (Okabe et al., 2009). However, iScience Article TROP2 high oval cells have been also reported to be poised toward developing into cholangiocytes (Seow et al., 2021). We found that EPCAM+ and TROP2+ populations were enriched for an adult biliary signature with a concurrent high expression of SOX9 (GSE102683; Figure 3A) (Ceulemans et al., 2017). This co-expression could indicate toward a possible enrichment of cholangiocytes in these samples, given that both TROP2 and EPCAM expression has been observed in cholangiocytes (Seow et al., 2021). Finally, progenitor cell samples isolated via side population techniques and progenitor cell niches isolated via laser microdissection showed relatively lower levels of c/EBPa, SOX9, adult hepatocyte and biliary signatures (GSE102683; Figure 3A) (Ceulemans et al., 2017). These cells instead showed an enrichment for TGFBR2 indicative of enrichment of progenitor cell populations in these samples (GSE102683; Figure 3A). Furthermore, when all these samples were analyzed together, the pairwise correlations between c/EBPa, SOX9, TGFBR2 and adult liver cell signatures showed consistent trends with what our model predicted (GSE102683; Figure S2B). Together, these observations suggest that the emergent dynamics of mechanistic model comprising c/EBPa, SOX9, and TGFb signaling can explain how these pathways act in concert to allow for different cell fates during various liver developmental and liver pathological scenarios.
Next, we assessed if gene expression changes observed in hepatic development and steatosis are seen in the context of liver reprogramming or during liver injury induced trans-differentiation. TGFb treatment can convert hepatocytes to cholangiocytes via trans-differentiation. Thus, we compared the TGFb signaling levels in a peripheral bile duct RNA-seq dataset from mouse livers (GSE108315) comprising cholangiocytes, hepatocytes, and hepatocyte-derived cholangiocytes. TGFb signaling was lower in primary hepatocytes, but enriched in primary cholangiocytes and those derived (trans-differentiated) from hepatocytes (Schaub et al., 2018). Consistent trends are noted for SOX9, while c/EBPa levels show the opposite trends, as expected ( Figure 3B). In another instance of trans-differentiation, in mouse oval cells (progenitor-like cells during liver regeneration) either derived from cholangiocytes or hepatocytes through chronic injury (Tarlow et al., 2014), the levels of SOX9 regulon and TGFb signaling were higher than in primary hepatocytes, but the adult hepatocyte signature was lower (GSE55552; Figure 3C). Finally, during reprogramming from human embryonic stem cells (hESCs) and/or hepatic stem cell line (HepaRG) (Dianat et al., 2014), hepatocytes -both mature or immature -exhibit low levels of TGFb signaling activity (compare fourth column with the sixth and seventh columns in Figure 3D), whereas cholangiocyte had higher or comparable levels of TGFb signaling activity as compared to hepatoblast cell state (compare second column with third column in Fig Previously, miR-337 has been implicated to promote a cholangiocyte fate during liver development (Demarez et al., 2018). Thus, we examined the differential gene expression programs induced by miR-337 mimic in in vitro cultures of immortalized hepatoblasts from E12.5 wild-type mice liver. c/EBPa regulon activity decreased significantly, but SOX9 regulon and TGFb signaling increased (GSE99890; Figure 3E), indicating the functional role of these players in committing to a cholangiocyte cell-fate. As another case study, we investigated hepatoblasts expressing adult stem cell marker LGR5 that have been shown to be capable of establishing both the hepatocyte and cholangiocyte colonies (Prior et al., 2019). We found a bimodal distribution of SOX9 regulon and TGFb signaling activity in the entire hepatoblast population as the developmental time progressed; however, no such shift in the levels was observed specifically within LGR5+ hepatoblasts ( Figure S2D), suggesting that a switch in TGFb signaling and SOX9 levels are required to commit the hepatoblast cells toward a cholangiocyte or a hepatocyte cell fate. Therefore, the analysis of these single cell and bulk RNA-seq datasets highlight how the emergent dynamics of c/EBPa, TGFBR2 and SOX9 as predicted by a minimal model signaling is capable of recapitulating the trajectories of cell-fate decisions as observed experimentally during hepatic development, regeneration and liver cell reprogramming.
Stochastic state switching and cellular reprogramming strategies among liver cell phenotypes based on gene regulatory network dynamics  (Sahoo et al., 2020a). We conducted stochastic simulations for the core gene regulatory network ( Figure 1A) (Kohar and Lu, 2018) and observed switching among cell fates for various parametric combinations. We observed that hepatocytes (high levels of c/EBPa) and cholangiocytes (high levels of SOX9) could switch between one another ( Figure 4A- iScience Article trans-differentiation as observed experimentally (Deng et al., 2018;Schaub et al., 2018). Similarly, we also observed that hepatoblasts (characterized by high levels of TGFBR2 and low levels of both SOX9 and c/EBPa) could switch to a cholangiocyte cell fate (characterized by high level of SOX9) and stably exist for relatively longer periods of time ( Figure 4A-middle panel). Finally, we also observed certain parameter sets showed dynamic transitions between all three cell states -hepatocytes, cholangiocytes and hepatoblasts ( Figure 4A-bottom panel), transitions which could not be seen in absence of noise ( Figure S4). These results indicate that multistable features of this gene regulatory network can confer cellular plasticity and consequent cell-fate switching in liver, which may be of paramount importance during injury repair and reprogramming.
After elucidating the dynamics of a core regulatory network for cellular decision-making during liver development and homeostasis, we predicted the cellular reprogramming efficiencies of different perturbations on this network and the final proportions of cell fates obtained. We performed either overexpression (OE) or down-regulation (DE) of c/EBPa, SOX9 and TGFBR2 individually and calculated the ratio of hepatocytes, cholangiocytes, progenitor-like cells and SOX9 + hepatocytes (see STAR methods for how the cell types have been defined for this analysis). Overall, we found that hepatocytes can be best enriched by either performing a c/EBPa OE or TGFBR2 DE ( Figure 4B). On the other hand, cholangiocytes can be best enriched for by SOX9 OE or TGFBR2 OE ( Figure 4C). For the progenitor-like state (high levels of TGFBR2; low levels of both c/EBPa and SOX9), one might expect that overexpression of TGFBR2 levels would lead to an enrichment of this state. However, we find out that TGFBR2 OE has a weaker impact in enriching progenitor cell population as compared to c/EBPa DE or SOX9 DE ( Figure 4D). This observation lends support to the hypothesis that to maintain a progenitor cell state population, both differentiation programs (cholangiocyte-inducing and hepatocyte-inducing) need to be inactive ( Figure 4E). This observation is likely to be more generic, given that various underlying gene regulatory networks for cell-fate decisions have similar topology as the one shown here (Duddu et al., 2020;Zhou and Huang, 2011). All above-mentioned trends remain qualitatively unchanged when TGFb was perturbed too. Further, SOX9 + hepatocytes, as expected, iScience Article were most enriched for by SOX9 OE instead of by TGFBR2 OE, although SOX9 is activated by TGFb signaling (Figures S3A and S3B), implying that enriching SOX9 + hepatocytes may need a more specific agonist facilitating SOX9 transcription, instead of a generic activation of TGFb signaling pathway. These simulations offer testable predictions for executing effective reprogramming strategies when liver cells of one phenotype switch to another sibling fate.
Overall, these observations endorse that concerted changes in levels of these 'master regulators' can orchestrate corresponding phenotypic plasticity and heterogeneity, as seen in vitro and in vivo.
Bifurcation analysis reveals spatial pattern formation and other dynamical properties associated with gene regulatory network in liver cells After investigating different possible phenotypes enabled by the network comprising c/EBPa, SOX9 and TGFBR2, we investigated the role of TGFb in controlling the transitions between hepatocyte and cholangiocyte cell states. For a set of kinetic parameters estimated from literature, we mapped different cellstates observed at various levels of TGFb through a bifurcation (dose-response curve) diagram (see STAR methods). At lower TGFb levels, we found c/EBPa can exist at two levels -high and low -mapping onto a hepatocyte and a progenitor/cholangiocyte fate respectively. Correspondingly, when c/EBPa levels were high, TGFBR2 and SOX9 levels were both low. On the other hand, when TGFBR2 and SOX9 were high, the levels of c/EBPa were lower, thus indicating bistability ( Figure 5A). As levels of TGFb increase, the levels of both TGFBR2 and SOX9 increase with a concurrent disappearance of the hepatocyte cell state (Figure 5A), which can be interpreted as the maturation of progenitor cells into a distinct cholangiocyte-like cell state and subsequent destabilization of the hepatocyte cell fate at higher levels of TGFb. In developmental contexts, it has been observed that the levels of TGFBR2 increase along with those of SOX9 as cholangiocytes develop/mature from a hepatoblast stage (Takayama et al., 2014), thus corroborating the bifurcation diagram seen here. This trend is also consistent with a graded mode in which cholangiocytes develop, i.e., not an abrupt change to a mature state (Yang et al., 2017). Furthermore, from the bifurcation diagram, we see that the progenitor cell state has higher level of SOX9 compared to hepatocytes, but lower than cholangiocytes ( Figure 5A). This observation can explain in vivo scenarios where progenitors having high levels of SOX9 can give rise to both biliary cells and hepatocytes in the context of liver development (Furuyama et al., 2011).
After examining the temporal dynamics of cell-fate commitment in the liver, we focused on spatial patterning of these phenotypes. Cholangiocytes are more abundantly found near the portal vein. On the contrary, hepatocytes are found much more abundantly on the parenchymal side (Clotman et al., 2005). Hepatic progenitor cells (similar to hepatoblasts/oval cells in our model) have been proposed to be sandwiched between the cholangiocytes and hepatocytes in the adult liver (Ko et al., 2020). Hepatoblasts in periportal region are likely to adopt a cholangiocyte fate and eventually form the bile ducts, whereas hepatoblasts in other regions of the lobule form hepatocytes (Raynaud et al., 2011). Such spatial segregation has been proposed to be a result of TGFb gradient (Ayabe et al., 2018;Clotman et al., 2005) and/or by mechanical cues induced Notch signaling variations in the periportal region (Kaylan et al., 2018). We probed whether we could reproduce this spatial patterning of hepatocytes/hepatoblast-like cells/cholangiocytes by simulations capturing a simple gradient of TGFb influence the gene regulatory network. Our simulations could reproduce the experimentally observed spatial pattern of phenotypes even without considering the local interactions between the neighboring cells in the simulation (i.e., the dynamics of each cell is independent of the neighboring cells). We demonstrated that under high levels of TGFb, a cholangiocyte cell state are more prevalent, whereas at intermediate and low levels of TGFb, hepatocytes and progenitor-like cells are seen to be abundant with the hepatocyte abundance increasing as the levels of TGFb drop ( Figure 5B). The analysis shows that spatial pattering is emergent from the gene regulatory network itself without necessarily invoking additional spatial interactions between neighbors. Our simulation observations are further corroborated by analysis of spatial gene expression profiles in liver tissue (GSE84498) (Halpern et al., 2017). We observed that the expression levels of c/EBPa drop toward the portal vein region with a concurrent increase in the levels of the SOX9 ( Figure 5C). These observations provide insights into the spatial organization of distinct cell types seen in liver development and homeostasis as a function of spatiotemporal dynamics of underlying gene regulatory networks.
A hallmark feature of multistable regulatory networks is hysteresis. To characterize the hysteretic behavior for this network, we first considered a ''pure population'' of hepatocytes and monitored, as a function of ll OPEN ACCESS iScience 25, 104955, September 16, 2022 9 iScience Article TGFb levels, the proportion of cells that switched to the progenitor/cholangiocyte state and the proportion of cells that did not switch states. Similarly, we started from a population of ''pure progenitors'' and tracked the relative proportion of cells in the three states (hepatocytes, cholangiocytes, and progenitor cells) as a function of TGFb levels. For a non-hysteretic system in nature, the proportions of cells in each of these three states would have been the same from the two starting populations (only hepatocytes or only progenitor cells) for a given level of TGFb. However, we observed distinctly different profiles for the different starting populations ( Figure 5D). For instance, at TGFb = 0, and starting from a population of only hepatocytes, the cells show minimal switching, instead maintain hepatocyte fate ( Figure 5D-left panel). But if one starts from a population of only progenitor cells, at TGFb = 0, only $30% switch to a hepatocyte fate and the remaining 70% retain their progenitor cell fate. Similarly, at TGFb = 15,000 molecules, if we start from a population of only progenitor cells, almost none of them switch to a hepatocyte fate but if we start from a population of only hepatocytes $20% of the cells will switch to a progenitor cell fate ( Figure 5D-middle panel). We also observed that the proportion of cholangiocytes obtained (when starting from a pure population of either hepatocytes or oval cells) has little effect on its starting population ( Figure 5D-right panel), which is expected, given that the cholangiocyte state is stable mostly beyond the bistable region where the hysteretic pattern is less pronounced. Although hysteresis in liver cell phenotypes requires experimental validation, it iScience Article has been experimentally witnessed for other multistable biological systems such as epithelial-mesenchymal plasticity and lactose metabolism (Celià -Terrassa et al., 2018;Ozbudak et al., 2004).
Finally, we also investigated possible trajectories of cellular reprogramming enabled by modulating the strength of individual edges in the network. We hypothesize that such alteration of strength of induvial edges can, in principle, happen via epigenetic regulations, but there is no current evidence for the same. Nonetheless, it is important to understand how such alterations might affect the dynamics of the gene regulatory network in question that can be later investigated experimentally. Given that TGFBR2 DE was shown to enrich for hepatocytes, we probed the impact of a stronger inhibition of TGFBR2 by SOX9 ( Figure 6A). For the control case, the bifurcation diagram showed that at TGFb =0, the system is bistable ( Figure 6A; blue curve in Figure 6B), i.e., in the absence of TGFb signaling, a certain fraction of progenitor-like cells will be present in the population, thus blocking a transition to all-hepatocyte population, when starting from a population of only cholangiocytes. However, on increasing the strength of suppression of TGFBR2 by SOX9, bifurcation diagram changed, now enabling a scenario of only hepatocyte cells at TGFb =0, with modest changes in steady state expression values of c/EBPa and SOX9 (green curve in Figure 6B). Further, cholangiocytes and hepatocytes can co-exist and likely stochastically switch even at higher values of TGFb signaling, without necessarily requiring transitioning via a hepatoblast state. Thus, a stronger inhibition of TGFBR2 by SOX9 provides a permissive situation for experimentally observed direct transdifferentiation among the two mature liver cell types (  iScience Article Next, we probed the impact of altering the link strength from SOX9 to c/EBPa. When this inhibition is weaker (case 1), the progenitor-like and cholangiocyte states are not observed irrespective of the value of TGFb (green curve in Figure 6D). In contrast, in case of a stronger inhibition (case 2), hepatocyte state disappears (orange curve in Figure 6D). Thus, by changing the link strength from SOX9 to c/EBPa, landscape of cellular plasticity can be transformed such that only one cell state ('attractor') exists ( Figure 6E), compromising on the ability to switch cell-fates in the liver. These observations suggest potent mechanisms for cells to lock into specific cell-fates as well as pinpoint how a delicate balance between different links in underlying regulatory network allow multistability and consequent phenotypic plasticity -a fundamental tenet of development and regeneration across organs and organisms.

DISCUSSION
Phenotypic plasticity is a ubiquitous phenomenon via which cells can exhibit diverse phenotypic states which can interconvert among themselves. It is thought to be a 'bet-hedging' strategy that can be implemented by cell populations to survive through various bottlenecks, such as drug-tolerant persisters seen ranging from microbial populations to cancer cells (Sahoo et al., 2020b;van Boxtel et al., 2017). Understanding the first principles behind cell-fate decision making in developmental systems can be paramount to elucidating mechanistic underpinnings of phenotypic plasticity in biological systems in general. Furthermore, such an understanding can propel rational strategies for reprogramming of cell types in vitro (Qian et al., 2018). The hepatocyte-cholangiocyte decision-making offers an ideal system, given their common progenitor (hepatoblasts) that can differentiate into these divergent cell-fate trajectories. Furthermore, once differentiated, these cell still retain the capacity to trans-differentiate to one another during liver injury (Gadd et al., 2020).
Here, we identified and analyzed a minimalistic gene regulatory circuit that can, in principle, explain various cellular phenotypes observed in the liver and their cell-fate decision under the influence of specific signaling cues. This network comprising SOX9, TGFBR2 and c/EBPa is multistable in nature; enabling multiple co-existing phenotypes that can switch in presence of biological noise. The experimentally observed spatial patterning of these cell types in the liver can also be explained as an emergent property of this gene regulatory network without explicitly considering spatial interactions between the neighboring cells in our minimalistic model. Finally, we explore possible perturbations in the gene regulatory network that can enrich for various cell phenotypes, driving reprogramming.
The in-silico model constructed here, like all model systems, has limitations. First, the hepatocyte and cholangiocyte phenotypes are proxied by individual 'master regulator' transcription factors, which makes it tricky to identify various heterogeneous 'micro-states' such as Axin+ hepatocytes, Tert High hepatocytes, hybrid periportal hepatocytes and Lgr5+ hepatoblasts (Li et al., 2020;Prior et al., 2019). Second, we have excluded several potentially important key regulators. For instance, we have not considered SOX4 as a distinct node in the gene regulatory network as its expression of SOX4 is largely correlated with that of SOX9 in multiple datasets we analyzed ( Figure S5). Including SOX4 in subsequent models can help identify 'micro-states' during development and disease. Similarly, we excluded c/EBPb, whose overexpression has been shown to promote cholangiocyte differentiation by enhancing TGFBR2 expression (Takayama et al., 2014). This observation, along with antagonistic trends in expression levels and/or regulation of c/EBPb and c/EBPa (Jakobsen et al., 2013;Takayama et al., 2014), suggests that c/EBPb can also be included in subsequent models to elucidate its role in development of cholangiocytes. Expanding the regulatory network to include these other nodes may alleviate this limitation; for instance, on including TGFb, our simulations showed the emergence of SOX9 + hepatocytes. These cells are reported to behave as bipotent progenitors after liver injury (Han et al., 2019;Tanimizu et al., 2017). Third, we only investigate the transcriptional control of these cell-fates; they can be affected by other layers of regulation such as epigenetics (Aloia, 2021). Future efforts to interrogate the impact of chromatin changes on cell-fate decisions (Jia et al., 2019) will be important. Fourth, our model cannot rigorously distinguish biological differences that might exist between developmental bipotent progenitors (hepatoblasts) and adult bipotent progenitors (oval cells) or between these progenitors and cells undergoing hepatocyte-to-cholangiocyte conversion (intermediate, and possibly metastable, states). Within the scope of this study, these diverse biological entities are considered to be similar, but they may be different in in vivo conditions. For instance, excessive activation of TGFb signaling in the liver is associated with appearance of ''hepatobiliary'' features; i.e. ''hybrid'' cells co-expressing both hepatocyte and cholangiocyte markers (Raynaud et al., 2011). Our analysis opens up a wide range of questions for which experiments can be designed to test out novel hypothesis. Based on the bifurcation analysis presented here and earlier literature, TGFb treatment is known to promote SOX9 levels by increasing activity of Notch/TGFb signaling (Wang et al., 2018), but it still remains to be seen if the prolonged treatment of TGFb antagonists can push the mature/immature cholangiocytes into a progenitor-like/hepatocyte cell fate. Furthermore, even if the program is reversible in the initial stages, is there a point of no return that prevents cholangiocytes from reverting back to hepatocytes/progenitor-like states, similar to those observed in other instances of phenotypic plasticity? (Tripathi et al., 2020) If yes, what mechanisms underlie a distinction between cell decision-making (reversible) and cell-fate 'locking' (irreversible). Juxtaposing stability versus plasticity (switching ability) in a multistable system can offer new insights into the operating principles invoked during canalization as well as organ injury repair. TGFb is a well-known inducer of a mesenchymal state in liver development and regeneration (Choi and Diehl, 2009;Karkampouna et al., 2012) as well as in cancer. However, cholangiocytes are known to be epithelial, even though TGFb appears to be crucial in their development. How the cholangiocytes maintain an ''epithelial phenotype'' under the influence of TGFb signaling remains to be answered. Potentially, interplay of TGFb with other signaling pathways such as Notch signaling can help resolve these context-specific discrepancies.
Overall, we demonstrate how a minimalistic gene regulatory network can provide insights in elucidating the dynamical properties of cell-fate decisions during liver development and regeneration. Our approach highlights common ''first principles'' that appear to be followed in varied biological contexts corresponding to healthy and diseased liver. Finally, our study provides a conceptual framework to guide cellular reprogramming strategies in pre-clinical and/or clinical setups by a better predictive understanding of emergent dynamics of networks driving liver cell-fate decisions.

Limitations of the study
Our minimalistic model explains how different cell types in the liver can switch among one another, but does not have the granularity to distinguish among certain biological micro-states or phenotypes seen in liver during development, or injury repair. Including other relevant master regulators to the network considered here would be critical in explaining such phenotypes. Also, other layers of biological regulation such as epigenetic, metabolic, and translational are not yet integrated into our framework to describe their role in liver cell plasticity. Finally, some of our model predictions and hypotheses would require further experimental validation in a context-specific manner.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:  iScience Article iScience Article depends on production rate of the node, the effect of regulatory links incident on the node (modeled as a shifted Hill's function (Lu et al., 2013)) and the degradation rate of the node. Each parameter in the set of unknown parameters for ODEs is randomly sampled from a biologically relevant range. The sampling of parameters is done so as to ensure that it generates a representative ensemble of models for a specific circuit topology. The range for production rate varies from 1 to 100 whereas the range for the degradation rates varies from 0.1 to 1 (arbitrary units). The fold change parameter associated with each link is assumed to be in the range of 1-to 100-fold. The Hills coefficients sampled are assumed to vary from 1 to 6. Furthermore, the threshold parameter is carefully chosen based on the gene network topology features such as inward regulations to ensure that each link in the network has a roughly equal chance of being functional or not functional. This ensures that steady state solutions obtained sample most if not all of the possible states allowed by the gene regulatory network without biasing the system towards a particular set of solutions.
After such sampling, the set of parameterized ODEs is solved to get different possible steady state solutions. The set of ODEs can be multistable, i.e., multiple sets of steady state concentrations satisfy the set of ODEs. The program samples 10000 different sets of parameters. For each parameter set, RACIPE chooses a random set of initial conditions (n = 100) for each node in the network and solves, using Euler's method, with the set of coupled ODEs that represent interactions among the different nodes in a network. For each given parameter set and initial conditions, RACIPE reports the steady-state values for each of the nodes in the network. The steady state values were then Z-normalized where the z-normalized expression value (z i ) is given by the term: where E mean and E std is the mean and standard deviation of the expression levels of a given node across all its steady state solutions.
The above procedure was followed to generate RACIPE results for the network topology given in Figure 1. The Z-normalized steady state values were then plotted as a clustermap as in Figure 1B. PCA was performed with the major clusters being labeled via hierarchical clustering setting n = 3 as in Figure 1C. Kernel density maps overlaid on top of histograms for steady state expression for each node was plotted as in Figure 1H. Scatterplots for steady state solutions for each pair of nodes was plotted as presented in Figure S1B. The variant circuit as shown in Figure S1C was also simulated and processed using the methods described above to generate Figures S1D and S1E.
The perturbation analysis in Figures 4B-4D and Figure S3A were done by performing RACIPE analysis on the system by either over expressing (OE) or down expressing (DE) the given node by 10-fold. The Z score normalization of these perturbation data was done with respect to the control case. Similarly, perturbation analysis was done on variant circuit shown in Figure S1C (shown in Figure S3B).
For spatial simulations, we started the system from an all-hepatocyte state and simulated 15 cells at each level of TGFb (0 to 40000 in steps of 2500). Same classification of biological phenotypes as above was used to assign cell types on the spatial pattern.

Data analysis
All bulk microarray and RNA-Seq pre-processed datasets were obtained from publicly available GEO datasets. Gene signatures for adult hepatocyte program, adult biliary program and adult hepatobiliary program were obtained from (Segal et al., 2019). The gene sets specifying c/EBPa regulon and SOX9 regulon were obtained from (Møller and Natarajan, 2020). Hallmark TGFb signaling pathway was obtained from MSigDB (Liberzon et al., 2011). Gene set activity were estimated via ssGSEA for all bulk samples and via AUCell for all single cell datasets (Aibar et al., 2017;Subramanian et al., 2005). For trajectory analysis the coordinates of each cell was obtained by performing a 3D PCA on the gene sets provided in (Segal et al., 2019).
We have used Shannon entropy as a proxy for the transcriptional diversity of the single cells in the dataset GSE90047. We calculate the Shannon entropy for the distribution of all genes' expression (in TPM) as a single quantity for a given cell.
For TGFb signaling, we included multiple metrics: expression levels of TGFB1, expression levels of TGFBR2, TGFb signaling activity ssGSEA scores. Representative cases for individual datasets are shown in corresponding main figures; all pairwise comparisons are given in Table S3. Most entries across datasets show the expected trend of TGFb signaling levels at a reduced level in the hepatocytes (fold-change levels <1 when compared to cholangiocytes and/or oval cells).

QUANTIFICATION AND STATISTICAL ANALYSIS
We computed the Spearman correlation coefficients and used corresponding pvalues to gauge the strength of correlations. For statistical comparison between groups, we used a two-tailed Student's t-test under the assumption of unequal variances and computed significance. Details of statistical analysis, definitions for significance and abbreviations can be found in figure legends.