NRF2-dependent Epigenetic Regulation can Promote the Hybrid Epithelial/Mesenchymal Phenotype

The epithelial-mesenchymal transition (EMT) is a cellular process critical for wound healing, cancer metastasis and embryonic development. Recent efforts have identified the role of hybrid epithelial/mesenchymal states, having both epithelial and mesehncymal traits, in enabling cancer metastasis and resistance to various therapies. Also, previous work has suggested that NRF2 can act as phenotypic stability factor to help stablize such hybrid states. Here, we incorporate a phenomenological epigenetic feedback effect into our previous computational model for EMT signaling. We show that this type of feedback can stabilize the hybrid state as compared to the fully mesenchymal phenotype if NRF2 can influence SNAIL at an epigenetic level, as this link makes transitions out of hybrid state more difficult. However, epigenetic regulation on other NRF2-related links do not significantly change the EMT dynamics. Finally, we considered possible cell division effects in our epigenetic regulation model, and our results indicate that the degree of epigenetic inheritance does not appear to be a critical factor for the hybrid E/M state stabilizing behavior of NRF2.


INTRODUCTION
Epithelial-Mesenchymal Transition (EMT) is often essential for various physiological and pathological processes such as wound healing, cancer metastasis, and embryonic development (Jolly, 2015). During EMT, cells tend to lose cell-cell adhesion and gain migration and invasion capabilities. Initially, most research assumed that EMT acted as a binary process, i.e., cells typically underwent full EMT, leading to independently-moving spindle-shaped cells. However recently, it has become clear that cells can also become stabilized in one or more hybrid epithelial/mesenchymal (E/ M) states, states which exhibit a combination of epithelial and mesenchymal traits and can lead to collective cell migration and enhanced tumorigenesis (Cheung and Ewald, 2016;Pastushenko and Blanpain, 2019;Tripathi et al., 2020a).
Nuclear factor E2-related factor 2(NFE2L2, i.e., NRF2) is a transcription factor which has been shown to play an important role in preventing cells from completing full EMT during injury-induced wound healing (Riahi et al., 2014). Previous mathematical models and knockdown experiments have proven that NRF2 can act as a "phenotypic stability factor" (PSF) which enables cells to more readily access a hybrid E/M state (Bocci et al., 2019). In these mathematical models, NRF2 interacts with the core EMT circuit via three regulatory links: 1) NRF2 inhibits Snail (Zhou et al., 2016); 2) E-cadherin inhibits the nuclear accumulation and transcriptional activity of NRF2 (Kim et al., 2012); 3) miR-200 promotes NRF2 by inhibiting Keap1 (Eades et al., 2011). All these interactions involving NRF2 act at a transcriptional or post-transcriptional level, but do not directly invoke epigenetic mechanisms such as DNA methylation or histone modification. These additional mechanisms have been shown to govern the extent and reversibility of cell-state transitions among the epithelial, mesenchymal and hybrid E/M phenotypes (Dumont et al., 2008;Jia et al., 2019;Jia et al., 2020;Serresi et al., 2021). Thus, the goal of this work is to ask how the aforementioned role of NRF2 can be modulated via epigenetics.
We integrated a mathematical model of coupled NRF2-EMT dynamics with a previously established phenomelogical epigenetic regulation framework (Miyamoto et al., 2015), and studied how epigenetic regulation could affect the stability of the hybrid E/M state. Specifically, we determined the effects of adding epigeneticsbased regulatory terms individually in all the three abovementioned NRF2 related-links. These terms dynamically modulate regulatory thresholds based on corresponding transcriptional activity (see Methods section). We found that incorporating epigenetic feedback affecting the inhibition of SNAIL, in other words, reducing the threshold for inhibition of SNAIL as a dynamic function of NRF2 levels, can significantly stablize the epithelial and hybrid E/M states. This effect was validated by a population dynamics analysis, and was not present for the other two links involving NRF2. Finally, we investigated the role of cell division in our model (Tripathi et al., 2020b), and our results suggest that epigenetic fluctuations due to cell division do not appear to play an important role in altering this stabilizing behavior of NRF2. FIGURE 1 | Epigenetic feedback-mediated dynamics of NRF2 coupled to core EMT circuitry. (A) A regulatory network for EMT that consists of two microRNA-TF mutually inhibitory circuits: miR-34/SNAIL and miR-200/ZEB. Signal I represents external EMT-inducing signals such as HGF, NF-κ B, Wnt, TGF-β and/or HIF1α. NRF2 module is added to the core networks by three distinct links. Blue links represent the inhibition from Keap1 and E-cad on NRF2 respectively, and the red link represents the inhibition on SNAIL from NRF2.(B-D) Bifurcation diagrams of miR-200 levels for the network shown in Panel 1A, with I as the bifurcation parameter. Solid lines represent stable states, dashed lines represent unstable states. Black lines correspond to the circuit without epigenetic regulation, and blue lines correspond to a circuit including epigenetic feedback. In (B), the epigenetic feedback is on the inhibition of SNAIL by NRF2. In (C), the epigenetic feedback is on the inhibitory link from E-cadherin to NRF2. In (D), the epigenetic feedback is on the inhibtion of NRF2 by Keap1.

Epigenetic Feedback on Inhibition of Snail by NRF2 May Stablize Hybrid E/M State
In a previous study, we incorporated a phenomenological treatement of epigenetic regulation into a microRNA-based model for a core EMT circuit consist of two interconnected mutually inhibitory loops between miR-34/SNAIL and miR-200/ZEB (Lu et al., 2013). We also analyzed additional effects that occur when the additional transcription factor GRHL2 is linked to the core model (Jolly et al., 2016;Jia et al., 2020) Here, we focus on an extension which adds to the core EMT model an explicit NRF2 module ( Figure 1A). Based on experimental data, the model assumes that NRF2 inhibits Snail, while E-cad and Keap1 inhibit NRF2. E-cad and ZEB are mutually inhibitory, and miR-200 inhibits Keap1. In this model, therefore, there are three links related to NRF2, i.e., NRF2-SNAIL, NRF2-E-cad and NRF2-Keap1.
We added epigenetic feedback individually to each link, and varied its strength in order to determine the range of possible effects. In each case, we first quantified these effects by deriving a new steady-state bifurcation diagram displaying how the allowed EMT states vary with an external inducing signal. In one of the three resultant bifurcation diagrams ( Figure 1B), where epigenetic feedback regulates the inhibition acting from NRF2 on SNAIL, the range of hybrid E/M state existence is greatly increased (miR-200 > 15,000 molecules: epithelial state; miR-200 levels < 5,000 molecules: mesenchymal state; 5,000 < miR-200 < 15,000 molecules: hybrid E/M state). Importantly, this increase is due to its "rightward" extension, which means even if the external EMT-inducing signal is very strong, the cell is still able to maintain its hybrid state and not automatically switch to a fully mesenchymal phenotype. In contrast, when epigenetic feedback is added to either of the two links which inhibit NRF2, the bifurcation diagrams are relatively unchanged as compared to the case when there is no epigenetic feedback ( Figures 1C,D). Based on these results, we focus below on the inhibition link from NRF2 to SNAIL.

Dynamical and Population Analysis Indicate the Stabilization of Hybrid E/M State
Based on results from the bifurcation analysis, we expect that including epigenetic feedback on the inhibition of SNAIL by NRF2 will make a significant difference in dynamical simulations which include noise, thereby allowing cell-state transitions to spontaneously occur (Sahoo et al., 2021;Li et al., 2016). Here, we start with a cell in the epithelial state, and the initial signal is set to be at the median of the tristable region (i.e., the co-existence region of E, M and hybrid E/M states) for the control (no epigenetic feedback) case. Specifically, the tristable signal region only ranges from 65,000-74,000 molecules, so we choose the initial signal to be 69,000 molecules in order to make it possible for the cell to switch among all three states. In a dynamical simulation without any epigenetic feedback (Figure 2A), most of the time the cell stays in a mesenchymal state (high ZEB, low miR-200). Intsead, if we include the epigenetic effect ( Figure 2B), the cell spends increased amounts of time in either a hybrid state (medium ZEB, medium miR-200) or an epithelial state (low ZEB, high miR-200), thus reflecting a significant change in the mean residence time (MRT) of the three phenotypes. These results are reminiscent of changes seen in MRT for these phenotypes upon including other "phenotypic stability factors" similar to NRF2 (Biswas et al., 2019). Here, the time unit ζ 100 h (simulation time) is the time factor used in the epigenetic feedback term (Supplementary Section S2).
These trajectory results have direct implications for the population dynamics. To see this, we performed population simulations with and without epigenetic feedback. Here, each cell had their initial phenotype as epithelial, and the initial signal is again set equal to the median of the tristable region for the bifurcation diagram obtained without any epigenetic feedback. In case of no epigenetic feedback, the population reaches a steadystate distribution of 4.5% epithelial, 7.8% hybrid and 87.7% mesenchymal ( Figure 2C). In comparision, when epigenetic feedback is added on the inhibition of SNAIL by NRF2, the population distribution becomes 52.4% epithelial, 31.4% hybrid and 16.2% mesenchymal ( Figure 2D). This dramatic increase in the population percentage of cells in the epithelial or hybrid states is consistent with the individual cell dynamics (Figures 2A,B).
So far, we have investigated the stabilization of hybrid and epithelial states. Following that, we now study the competition between hybrid state and mesenchymal state. We chose the mean of signal to be 140K, which is mainly in the {E/M, M} region for models both with or without epigenetic feedback cases, and which lies close to the edge of the epithelial state in the bifurcation diagram ( Figure 3A). We implemented a population analysis for 100 cells. Without epigenetic feedback, the percentage of hybrid E/M cells continues decreasing when the percentage of mesenchymal is increasing, and the population reaches a distribution consisting of 5%E, 10% E/M, 85% M at the end of our simulation ( Figure 3B). However, with the epigenetic modification, the population maintains a high level of hybrid, and still has 84% E/M in the end, while the E and M stabilize around 8% respectively ( Figure 3C). Given the huge difference revealed by this analysis, we suggests that epigenetic feedback by NRF2 may help cells maintain a hybrid state, when mainly competing with mesenchymal fates.

Epigenetic Feedback by NRF2 Helps Cells Maintain Epithelial State When Competing With Hybrid State
Next, we investigated the competition between epithelial and hybrid states. In order to increase the range of the tristable region, we added a GRHL2 module to our baseline model ( Figures  4A,B). We then treat cells with a relatively high external signal (I 170 K to 150 K, green dotted line, {H,M} region for epigenetic model, {M} region for the non-epigenetic model) for different time durations, and subsequently reduce the signal to a lower level (I 102 K, purple dotted line, {E,H,M} region for both the epigenetic and non-epigenetic cases). When the signal is high for only a short time period in circuits with or without epigenetic feedback, the cells typically remain in their initial epithelial state ( Figure 4C(1), (3)), though cell may occasionally transit to hybrid state in cases without epigenetic feedback. The population results indicate a final stable distribution of 30% hybrid E/M, 70% E when without epigenetic feedback ( Figure 4C(2)), and fully epithelial states with epigenetic feedback. When the duration of the high signal is extended in the model without epigenetic influence, the cells can reach hybrid state and stay for a certain time period ( Figure 4D(1), (3)). In some samples with epigenetic feedback, cell still stay in the epithelial state ( Figure 4D(3)). Over a long time period, the population analysis shows that in both cases the percentage of hybrid state is slowly decreasing while the percentage of epithelial state is increasing ( Figure 4D(2), (4)). Meanwhile, the increase in percentage of cells in an epithelial state is faster in the case of an epigenetic feedback, and the whole population eventually has a higher percentage of epithelial cells, as compared with the no epigenetic feedback group. This comparison also indicates that the system requires a longer time to re-organize its population when epigenetic feedback plays a role. Finally, if we keep increasing the treatment duration of high signal in the baseline case, a cell can complete the full epithelial-mesenchymal transition, and maintain a mesenchymal state after signal reduction ( Figure 4E(1)); the whole population becomes fully mesenchymal in this case as the noise level is insufficient to cause transitions. However, the population can still maintain hybrid state cells when adding epigenetic feedback ( Figure 4E(3)), and the distribution is similar to that of Figure 4D(4). These results reveals that when epigenetic feedback regulates the inhibition from NRF2 to SNAIL, the epithelial state becomes more stable when it competes with hybrid states, and also the epigenetic feedback increases the re-organizing time needed for our circuit. We again see the relative suppression of more mesenchymal states in the prersence of epigenetic effects.

Epigenetic Inheritance may not be Precisely Required for the Effect of NRF2 on Hybrid E/M Stability
Finally, we combined a previous cell division model with this epigenetic regulation model (Miyamoto et al., 2015;Tripathi et al., 2020b), in order to see how imprecise epigenetic memory would affect the population behavior and influence the behavior of NRF2 in terms of stabilizing hybrid E/M phenotype. In our original cell division model, the signal I_ext acquires a noise term during cell division, due to unequal distribution of signaling molecules; thus, the daughter cells may have different cell type as compared with the parent cell. Here, we considered the model containing the aforementioned epigenetic term governing SNAIL's inhibition by NRF2, and introduced noise into the threshold value in the Hill function which governs the inhibition from NRF2 to SNAIL, during each cell division. This stochasticity mimics the possibility that epigenetic marks may not be perfectly reproduced in daughter cells. We simulated three cases with threshold noise values 0, 100 and 300 K (intial threshold 1000 K). Surprisingly, all cases gave rise to similar distributions, around 50% epithelial, 13% hybrid and 37% mesenchymal ( Figure 5A). In other words, the failure of precise epigenetic inheritance seems to have a negligible effect on the population dynamics of EMT.
In order to better understand this result, we plotted the bifurcation diagram for the threshold as a control parameter, because our simulations varied the threshold value ( Figure 5C). We picked the value of signal I 110 K that lies in the tristable region. Then, we ran dynamical simulations to capture trends in individual cells, keeping all intial parameters fixed at the epithelial state except for the threshold value. We found that only varying this threshold is not sufficient for cells to make transitions, even when the initial threshold is set to equal the value for the stable mesenchymal state ( Figure 5B). In fact, the epigenetic dynamics is driven by the values of other variables which quickly return to their value in the epithelial state. Based on this analysis, we conclude that in this model, noise in the epigenetic state is not important in terms of determing the population structure.

DISCUSSION
Many computational and experimental studies have focused on a set of phenotypic stability factors such as GRHL2, OVOL2 and Np63α, which can prevent cells from undergoing the full EMT process and thus help cells maintain a hybrid epithelial/ mesenchymal state (Roca et al., 2013;Watanabe et al., 2014;Dang et al., 2015;Hong et al., 2015;Jolly et al., 2017;Bui et al., 2020;Silveira et al., 2020). NRF2 is also one of these phenotypic stability factors, and previously constructed mechanism-based NRF2-included EMT model successfully showed that NRF2 can stablize hybrid states (Bocci et al., 2019).
Here, we add additional regulatory dynamics to incorporate epigenetic effects. Our epigenetics-enhanced EMT model predicts that when there is epigenetic regulation modulating the inhibition of SNAIL by NRF2, the stabilization of hybrid states can be further enhanced. In such a case, cells can favor epithelial or hybrid E/M states, depending on the bifurcation region they lie in. Considering the epigenetic regulation, the epithelial state is more  Figure 4C). (C) Bifurcation diagram for x0ms based on inpit signal I, for the model including both epigenetic regulation on SNAIL's inhibition from NRF2 and NRF2's inhibition from Keap1.
Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 828250 7 stable when competing with the hybrid state, while hybrid state is greatly favored when competing with the mesenchymal state.
We also studied the implications of our phenomelogical model of epigenetic regulation for cell division. In our simulation, the epigentic feedback regulates cells in similar timescale compared with cell division timescale. Our results indicate that error-free epigenetic inheritance may not be critical in determing EMT population balance. This does not necessarily imply that cell division is not necessary, as cell division also introduces noise in all TF and microRNA concentrations (Baptista and Ribeiro, 2020). In fact, recent experiments (den Hollander et al., 2021) seem to indicate that suppressing division does interfere with EMT, albeit without demonstrating a specific causal connection. We should note that our result may be model-dependent, and approaches in which epigenetic effects are more directly coupled to the allowed cell states (Zhao et al., 2021) might give different results. We plan to further investigate this issue in future.
Our work adds to the growing literature on mechanism-based models to understand the emergent dynamics of partial and complete EMT (Hari et al., 2020), which complements approaches based on data-based or statistical models inferring the trajectories of EMT in a high-dimensional landscape (Sha et al., 2020;Sha et al., 2021). The influence of cell-cell communication with other tumor and/or stromal cells on EMT and the impact of EMT on therapeutic ressistance remain another active research area of investigation (Li et al., 2019;Scott et al., 2019;Bocci et al., 2021).

METHODS
The dynamics of our EMT regulatory circuit were simulated by Ordinary Differential Equations (ODEs), and all the mathematical equations and corresponding parameters are shown in Supplementary Section S1 Supplementary Table  S1. The simulations were implemented by Matlab.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.