Coupled Feedback Loops Involving PAGE4, EMT and Notch Signaling Can Give Rise to Non-Genetic Heterogeneity in Prostate Cancer Cells

Abstract

Non-genetic heterogeneity is emerging as a crucial factor underlying therapy resistance in multiple cancers. However, the design principles of regulatory networks underlying non-genetic heterogeneity in cancer remain poorly understood. Here, we investigate the coupled dynamics of feedback loops involving (a) oscillations in androgen receptor (AR) signaling mediated through an intrinsically disordered protein PAGE4, (b) multistability in epithelial–mesenchymal transition (EMT), and (c) Notch–Delta–Jagged signaling mediated cell-cell communication, each of which can generate non-genetic heterogeneity through multistability and/or oscillations. Our results show how different coupling strengths between AR and EMT signaling can lead to monostability, bistability, or oscillations in the levels of AR, as well as propagation of oscillations to EMT dynamics. These results reveal the emergent dynamics of coupled oscillatory and multi-stable systems and unravel mechanisms by which non-genetic heterogeneity in AR levels can be generated, which can act as a barrier to most existing therapies for prostate cancer patients.

1. Introduction

Phenotypic plasticity, i.e., the ability of cells to switch back and forth reversibly among different states (phenotypes), is a universal feature of adaptation to varying environments encountered by various biological systems [1]. This theme has been investigated in developmental and evolutionary biology in detail [2,3], and is gaining importance in the context of disease progression as well [4,5,6,7]. Further, this theme is well-studied in cases of bacterial and yeast populations [8,9,10,11,12], and is increasingly being investigated for mammalian cells as well [13,14,15]. Biochemical mechanisms underlying phenotypic plasticity and consequent non-mutational or non-genetic heterogeneity, and their implications in determining the fitness of individual cells and entire cell populations remain to be comprehensively elucidated [16,17,18,19,20,21]. Recently, phenotypic plasticity has emerged as an important player in facilitating resistance against many standard chemotherapeutic drugs and targeted therapies for multiple cancers [22,23,24]. Similar to persisters in a bacterial population, drug-tolerant persisters have been observed in cancer [24]. Thus, besides the well-studied genetic/genomic heterogeneity in drug resistance, phenotypic (non-genetic) heterogeneity can enable adaptive cross-drug tolerance for cancer cells [25,26,27]. Unlike genomic changes that are “hard-wired” and can be inherited by cell division, phenotypic changes are reversible and stochastic and may or may not be transmitted to the next generation [28,29,30,31].

Stochasticity in phenotypic plasticity and/or heterogeneity is a direct consequence of limited number effects of molecules involved in various biochemical reactions, including transcription, translation and cell-cell communication signaling [32,33,34]. Stochasticity can also drive cell-state transitions in a multistable system; this phenomenon may drive tumor repopulation after a phenotypic subpopulation has been selectively killed by the drug(s) [35]. Another regulatory level at which stochasticity can lead to cell-to-cell heterogeneity is promiscuity in protein-protein interactions. Such “conformational noise” is manifested particularly by intrinsically disordered regions/proteins (IDRs/IDPs); many oncogenes/tumor suppressor genes have been shown to contain IDRs [36,37,38,39]. While both these modes of biological noise have been probed separately, their combined effect on cancer cell dynamics remains elusive. Here, we investigate the case of prostate cancer through the lens of these two mechanisms.

In prostate cancer (PCa), androgen-deprivation therapy (ADT) has been a standard-of-care treatment for over 75 years. Resistance to ADT eventually occurs in most patients, leading to metastatic castration-resistant prostate cancer (CRPC) [40]). Progression after ADT has been often connected to the epithelial–mesenchymal transition (EMT) [41,42], a program that acts as the fulcrum of phenotypic plasticity. During EMT, epithelial cells tend to lose cell-cell adhesion and apicobasal polarity while gaining motility typical of mesenchymal cells. However, EMT in a tumor micro-environment, is neither irreversible nor a binary process. Such multistability in EMT can be driven by mutually inhibitory feedback loops among transcription factors and microRNAs whose relative levels dictate the cellular phenotype [43]. Multistability can lead to transitions among multiple cell states—epithelial (E), mesenchymal (M) and hybrid E/M—as observed for prostate cancer cells [44]. On the other hand, intrinsic disorder in the cancer/testis antigen PAGE4, has been suggested to regulate the signaling of androgen receptor (AR), a target of ADT. Differently phosphorylated versions of PAGE4 can form a negative feedback loop involving AR, which has been predicted to give rise to oscillations [45,46], thereby generating non-genetic heterogeneity in the levels of AR in an isogenic population. The existence of multiple molecular programs that induce phenotypic plasticity and heterogeneity raises interesting and so far, unanswered questions about their mutual interconnections, emergent dynamics, and clinical implications in the context of PCa. AR can form a mutually inhibitory loop with ZEB1, an EMT-inducer [41], thus coupling a multistable system with an oscillatory one. To shed light onto the emergent dynamics of the PAGE4 and EMT signaling networks, here we simulate the coupled dynamics of the AR-EMT circuit and demonstrate that, depending on the strength of bidirectional coupling between ZEB1 and AR, multistability can be seen in AR signaling and/or some oscillations can be seen in EMT circuitry. Besides intracellular signaling, we also investigate the effect of this coupling on Notch signaling, a cell-cell communication pathway. Our results reveal how AR signaling can display different nonlinear emergent dynamics (oscillations, bistable) and therefore generate and maintain non-genetic heterogeneity in an isogenic cancer cell population.

2. Results

2.1. Oscillations and Multistability in PAGE4-AR and EMT Standalone Signaling Networks

We have previously investigated standalone dynamics of PAGE4-AR and EMT circuits [46,47]. The PAGE4-AR circuit consists of three relevant PAGE4 phospho-forms (WT-PAGE4, HIPK1-PAGE4, CLK2-PAGE4), together with c-Jun and androgen receptor (AR). HIPK1 and CLK2 are both enzymes that can phosphorylate PAGE4 at different residues. The potentiation of c-Jun by HIPK1-PAGE4 can eventually drive the hyper-phosphorylation of HIPK1-PAGE4 to a short-lived CLK2-PAGE4 (Figure 1(Ai)). These interactions lead to a delayed negative feedback loop that can drive oscillations in this circuit with a typical period of 1–2 weeks [46]. This hyperphosphorylation of HIPK1-PAGE4 to CLK2-PAGE4 involves AR which is inhibited by c-Jun potentiation and can inhibit CLK2 upregulation. To further generalize our previous model, we investigate the circuit’s behavior for a variable strength of this negative feedback based on a fold-change parameter (${\lambda }_{PAGE4}$) that can be continuously varied between 0 and 1. ${\lambda }_{PAGE4}=0$ indicates maximal inhibition on AR (by c-Jun) and CLK2 (by AR). On the other hand, ${\lambda }_{PAGE4}=1$ indicates no inhibition for either of the two links at all, thus, the feedback loop is broken (see Methods). While a strong negative feedback leads to sustained oscillations of AR (Figure 1(Aii)), a weak negative feedback leads to steady AR levels (Figure 1(Aiii)).

On the other hand, a core EMT circuit can be characterized by a mutual inhibition between the microRNA miR-200 and transcription factor ZEB. ZEB here denotes a family of transcription factors whose members are ZEB1 and ZEB2, and drive EMT, i.e., it is a EMT-inducing transcription factor, while miR-200 overexpression can drive the reverse of EMT—MET (mesenchymal–epithelial transition) [48]. ZEB can also self-activate directly/indirectly [49] (Figure 1(Bi)). Moreover, SNAIL acts as an external signal driving EMT by activating ZEB and inhibiting miR-200 [47]. Thus, a bifurcation diagram of miR-200 levels with respect to SNAIL displays transition from an epithelial (E: high miR-200, low ZEB) to a hybrid epithelial/mesenchymal (E/M: medium miR-200, medium ZEB) to a mesenchymal (M: low miR-200, high ZEB) state. Therefore, the standalone SNAIL/ZEB/miR-200 circuit can behave as a monostable, bistable, or tristable system with co-existence of E, E/M and M phenotypes depending on the level of SNAIL-driven EMT induction (Figure 1(Bii,iii)).

Next, in the following sections, we investigate how different couplings between the AR and EMT circuits can modify their standalone circuits’ dynamics and give rise to different states. First, to gain an understanding of how the AR circuit responds to external signals, we study a generic coupling between AR and an external node (X). Next, we explicitly couple the PAGE4 and EMT circuit through AR-ZEB mutual inhibition.

2.2. Negative Feedback Loops between AR and External Signals Give Rise to Oscillations, Monostability or Bistability

Before coupling the PAGE4-AR circuit with the EMT circuit, we first add a single node X to the PAGE4-AR circuit to understand how perturbations to AR signaling can modify its oscillatory dynamics (Figure 2A). We also allow X to self-regulate (either self-inhibit or self-activate), as seen for various “master regulators” of cell-state transitions [50,51]. This coupling between AR and X mimics the scenario of mutual inhibition between AR and ZEB1, and possible self-activation of ZEB1 [49]. We investigated the dynamics of this extended circuit at varied strengths of coupling between AR and X (represented by ${\lambda }_{DNFL}$), and different strengths of interaction in PAGE4-AR feedback loop (represented by ${\lambda }_{PAGE4}$). These values vary between 0 and 1, and the smaller the value, the stronger the inhibition (see Methods). In this circuit, ${\lambda }_{DNFL}$ represents the fold-change (and hence the strength of regulation) for both links: the inhibition of AR by X and vice versa.

At high ${\lambda }_{PAGE4}$ and low ${\lambda }_{DNFL}$, i.e., when AR and X inhibit each other strongly but the internal interactions in the PAGE4-AR circuit are weak, the system shows a bistable behavior—co-existence of (high AR, low X) and (low AR, high X) states, typical of double negative feedback loops [52,53] (Figure 2B; center panel, top left region). Conversely, at low ${\lambda }_{PAGE4}$ and high ${\lambda }_{DNFL}$, i.e., when AR and X inhibit each other moderately but PAGE4-AR circuit features a strong negative feedback loop, the system displays sustained oscillations, typical of delayed negative feedback loops [54,55,56] (Figure 2B; center panel, bottom right region). Interestingly, at both (high ${\lambda }_{PAGE4}$, high ${\lambda }_{DNFL}$) and (low ${\lambda }_{PAGE4}$, low ${\lambda }_{DNFL}$), AR is mono-stable, albeit at two different steady state levels (Figure 2B, center panel; top right and bottom left regions). In the former case, AR saturates at a higher level, perhaps due to lack of inhibition by X and/or other components of PAGE4-AR circuit. Conversely, in the latter case, AR saturates to a lower level. Thus, the coupled PAGE4-AR-X circuit can show bistable, monostable or oscillatory dynamics depending on relative strengths of the regulatory links.

Next, we study the effects of self-regulation on X, denoted by the fold-change parameter ${\lambda }_{XtoX}$. A value greater than 1 implies self-activation of X, whereas a value smaller than 1 implies self-inhibition of X. In the case of self-inhibition (${\lambda }_{XtoX}=0.1$), the bistable phase disappears and the system exhibits either monostable dynamics or oscillations (Figure 2C). This trend is consistent with observations that positive feedback loops facilitate multistability [43]. Oscillations are noted only at strong coupling within the PAGE4-AR feedback loop (i.e., smaller values of ${\lambda }_{PAGE4}$) (Figure 2C, yellow-shaded region); for weak PAGE4-AR coupling (i.e., higher values of ${\lambda }_{PAGE4}$), AR saturates at high steady state values, irrespective of the coupling strength between AR and X (Figure 2C, blue-shaded region). Conversely, in case of self-activation of X, the system largely behaves as in the case of no-self regulation, although with an increased parameter regime for bistability (Figure 2D). Thus, the self-regulation of X can alter the parameter regions which enable monostable, bistable and oscillatory dynamics for the PAGE4-AR-X circuit.

To further gain confidence in the abovementioned observations, we investigated the effect of altering other model parameters. All non-phosphorylation reactions here are modeled via shifted Hill function which describe the production fold-change of a given species as a function of inducer/inhibitor levels (see Methods); specifically, a Hill coefficient (n) quantifies how nonlinearly or steeply the fold-change depends on inducer/inhibitor level. We observed that the higher the value of n for AR-X coupling ($n_{DNFL}$) and the lower the value of n for coupling within PAGE4-AR circuit ($n_{PAGE4}$), the larger the parameter region enabling bistability. Conversely, lower values of $n_{DNFL}$ and/or higher values of $n_{PAGE4}$ drives oscillations (Figure S1A). These observations endorse our observations that the feedback loop of AR with X and that with PAGE4 and its phosphorylated forms push the system to behave differently: bistability in the former case, oscillations in the latter. Further, we have so far considered identical parameters for the shifted Hill functions denoting the inhibition of AR by c-Jun and that of CLK2 by AR. For the purpose of sensitivity analysis, even when we considered them to be independent parameters, we noticed that both these links have to be strong (i.e., λ$~0$) to facilitate oscillations in the levels of AR (Figure S1B). Put together, all these results suggest that a strong “external coupling” (i.e., double negative feedback loop between AR and X) favors bistability, while a strong “internal coupling” (i.e., negative feedback loop formed between AR and phospho-forms of PAGE4) drives oscillations in PAGE4-AR-X coupled circuit.

2.3. Dynamics of Coupled PAGE4-AR-EMT Circuits

After characterizing the PAGE4-AR response to generic external signals, we investigate the dynamics of coupled PAGE4-AR and EMT circuit. In this coupled circuit, AR and ZEB1 inhibit each other (similar to the generic AR-X coupling explored in the previous section), while SNAIL (S) is an external EMT-inducer (Figure 3A). To investigate how the AR-ZEB coupling modifies the dynamics of the coupled circuit, we probe the system’s dynamics for various strength combinations of the AR-to-ZEB inhibition (described by the fold-change parameter ${\lambda }_{AtoZ}$) and ZEB-to-AR inhibition (described by the fold-change parameter ${\lambda }_{ZtoA}$). Given that the coupled PAGE4-EMT circuit includes time delay and can potentially give rise to oscillations, we inspect the behavior of the circuit by evaluating its temporal dynamics for various combinations of coupling strengths and SNAIL level.

As a first step, we study the effects of weak bidirectional coupling, i.e., ${\lambda }_{AtoZ}$ = ${\lambda }_{ZtoA}$ = 0.9. As expected, here, both circuits largely show their standalone dynamics, i.e., oscillations for PAGE4-AR circuit and monostability/multistability for the EMT circuit depending on the levels of SNAIL. Interestingly, the EMT circuit can show oscillations of very small magnitude on top of its steady states obtained, especially at SNAIL levels enabling multistability (Figure S2).

Next, we consider the case where ZEB inhibits AR strongly while AR only weakly inhibits ZEB (${\lambda }_{ZtoA}=0.1$, ${\lambda }_{AtoZ}=0.9$). In this regime, ZEB can be approximated to behave as an external input to AR without a strong feedback. For low SNAIL values (S = 160 K), AR oscillates, whereas miR-200 relaxes to a high value typical of the epithelial state (Figure S3A). This observation can be attributed to relatively weaker induction or overall low levels of ZEB which are insufficient to reduce the levels of either miR-200 or AR. As SNAIL levels are increased (S = 195 K), the circuit can attain two possible stable states—an epithelial (high miR-200) state with AR oscillations or a mesenchymal (low miR-200) state with low AR (Figure 3B). On further increasing SNAIL levels (S = 200 K), a third hybrid E/M (medium miR-200) state emerges. (Figure 3C). While the epithelial state allows AR oscillations, a partial or complete EMT (i.e., hybrid E/M or M state) tends to have suppressed or no oscillations. This difference can emerge due to varied levels of ZEB in these states; higher ZEB levels can inhibit AR and hence dampen the oscillations seen in PAGE4-AR coupled circuit. Further increase of SNAIL levels (S = 215 K) eliminates the epithelial state, and the EMT circuitry shows coexistence of the hybrid E/M and mesenchymal states (Figure S3B). An even higher value of SNAIL (S = 240 K) drives a monostable mesenchymal phase in the EMT bifurcation diagram (Figure 1B). For these values of SNAIL, oscillations in AR disappear and the AR dynamics converge to low steady state values, given the higher levels of ZEB and consequently a strong inhibition of AR by ZEB (Figure 3C, D). Overall, in the case of strong inhibition of AR by ZEB, the multistability of the EMT circuit can be propagated to PAGE4-AR circuit, and SNAIL-driven EMT induction can dampen the PAGE4-AR oscillations.

The above-mentioned results consider a case of strong PAGE4 ‘internal coupling’, i.e., in the parameter region where the standalone dynamics of PAGE4-AR circuit is oscillatory $({\lambda }_{PAGE4}=0.1)$. Conversely, in the case of weak internal coupling $({\lambda }_{PAGE4}=0.9)$, AR saturates to a high steady state value (monostable) at low SNAIL values (S = 180 K) (Figure 3E). This difference can be explained by less ‘resistance’ offered by PAGE4-AR circuit in showing bistable behavior which can be overcome by relatively lower levels of SNAIL or ZEB. Increasing SNAIL levels (S = 190 K), however, induces bistability in AR levels, such that depending on the initial condition, cells can converge to either an (epithelial, high AR) state or a (mesenchymal, low AR) state (Figure 3F). Thus, multistability is passed on from EMT circuit to the PAGE4-AR circuit in this case too. Interestingly, the weak oscillations seen for miR-200 in the epithelial state for strong “internal coupling” (Figure 3C) also disappear in the case of weak “internal coupling”, because the PAGE4-AR circuit does not exhibit sustained oscillations in this case.

Next, we examine the case when we interchange the interaction strengths of both arms of the feedback loop between AR and ZEB1, i.e., when AR inhibits ZEB strongly, but ZEB inhibits AR weakly (${\lambda }_{AtoZ}=0.1$, ${\lambda }_{ZtoA}=0.9$). In this regime, AR is similar to an external signal to EMT circuit. Thus, it oscillates for any value of SNAIL due to the weak effect of ZEB on AR (Figure 4(Ai); Figure S3C). Conversely, the dynamics of the EMT circuit is largely driven by AR. Interestingly, miR-200 does not show multistability but rather oscillates with an amplitude that increases as a function of SNAIL (Figure 4(Aii–iv)).

Finally, we consider the case of strong coupling on both arms, i.e., when AR and ZEB1 both inhibit each other strongly (${\lambda }_{AtoZ}$ = ${\lambda }_{ZtoA}$ = 0.1). In this regime, both circuits strongly influence each other’s dynamics. For low SNAIL values (S = 160 K), AR oscillates and miR-200 assumes an epithelial level (Figure 4B). However, at higher SNAIL values (S = 200 K), two states become accessible: an epithelial state (high miR-200) with AR oscillations and a mesenchymal state (low miR-200) with low and steady AR levels (Figure 4C). Interestingly, AR oscillations propagate to miR-200 in the epithelial state. A further increase of SNAIL (S = 215 K) introduces a third, hybrid E/M state characterized by medium miR-200 levels and AR oscillations (Figure 4D). Finally, at very high SNAIL levels (S = 240 K), only a low miR-200, mesenchymal state is accessible, and correspondingly AR saturates to a low steady state value (Figure 4E). Overall, when the bidirectional signaling between ZEB and AR is strong, epithelial and hybrid E/M states are susceptible to AR oscillations, whereas a completely mesenchymal state suppresses AR plasticity.

To contextualize our results with respect to the model’s parameters, we further perform sensitivity analysis by vary each of the model parameters (one at a time) by +10% or −10%. We quantify parameter sensitivity in terms of variation for amplitude and period of AR oscillations (Figure S4). Generally speaking, oscillations are robust to local parameter change for both the cases weak effect of Zeb on AR (${\lambda }_{AtoZ}=0.1$, ${\lambda }_{ZtoA}=0.9$) i.e., Figure 4(Ai) and the strong effect of Zeb on AR (${\lambda }_{AtoZ}$ = ${\lambda }_{ZtoA}$ = 0.1) (Figure 4(Bi)). Parameters of production and degradation rate of AR and CLK2 are the most sensitive ones, which can be attributed to their effect in modulating the overall strength of the delayed negative feedback loop that generates oscillations in AR.

To summarize the results, we drew phase plots showcasing the accessible EMT states (Figure 5(Ai–iii)) and amplitude of AR oscillations (Figure 5(Bi–iii)) as a function of strengths of AR-to-ZEB and ZEB-to-AR inhibition (${\lambda }_{AtoZ}$ and ${\lambda }_{ZtoA}$), for three different values of SNAIL. By comparing the regions of the stability of EMT phenotypes and AR oscillation amplitude, it is possible to identify some general trends. First, monostability of the epithelial state corresponds to maximal amplitude of AR oscillations (blue region in Figure 5(Bi) overlaps with red region in Figure 5(Ai)). Conversely, AR does not typically oscillate with a large amplitude in regions where the mesenchymal state is the only state present, due to high levels of ZEB and consequently low levels of AR (red region in Figure 5(Biii) overlaps with slate region in Figure 5(Aiii)). The transition between oscillating and non-oscillating regimes seems fairly continuous, and intermediate oscillation amplitudes are found in multistable regions where various combinations of E, hybrid E/M and M states are accessible.

2.4. Integrating Notch–Delta–Jagged Signaling Circuit with Combined PAGE4-EMT Circuit

So far, we have considered phenotypic heterogeneity emerging from intracellular signaling. Communication between cells, however, can potentially modulate the PAGE4-EMT dynamics and propagate heterogeneity to other cells. Specifically, Notch signaling is widely recognized as a crucial mediator of cancer progression that operates via ligand-receptor binding between neighboring cells [50]. Notch signaling is activated by transactivation of Notch receptor by Delta and/or Jagged ligands presented by neighboring cells. This activation leads to cleavage of Notch, thus forming Notch Intra-Cellular Domain (NICD) which can translocate to the nucleus and affect the expression of multiple target genes, including Notch, Delta and Jagged themselves. NICD activates Notch and Jagged, but inhibits Delta (Figure 5A, Notch circuit). Therefore, Notch–Delta–Jagged signaling can lead to different kinds of feedback loops across cells. Notch–Delta signaling typically behaves as an intercellular toggle switch leading to (high Notch, low Delta) and (low Notch, high Delta) states–often called “lateral inhibition”. The (high Notch, low Delta) state is referred to as “Receiver” (R) because Notch is a transmembrane receptor and Delta is its ligand; similarly, the (low Notch, high Delta) state is referred to as a “Sender” (S). Notch-Jagged signaling, on the other hand, can lead to “lateral induction”, i.e., both neighboring cells exhibit (high Notch, high Jagged) levels and can exhibit a hybrid sender/receiver (S/R) phenotype. These different feedback loops mediated by Notch affects cell patterning and consequent cell-fate determination across biological contexts [50].

Previously, we have studied the coupled dynamics of EMT and Notch that were capable of exhibiting a maximum of four distinct states [57]. Drawing from these studies, here, we combined the three networks (PAGE4-AR, EMT, Notch–Delta–Jagged) at a single-cell level. Notch signaling is activated by transactivation of Notch receptor by Delta and/or Jagged ligands, thus we investigate the dynamical behavior of the coupled circuit at varied levels of such ligands. In other words, we model the coupled PAGE4-EMT-Notch circuit in an individual cell that is exposed to a varying external levels of Notch receptors (Notch) and ligands (Delta, Jagged). These external receptors and ligands mimic the effect of neighboring cells and activate the intracellular Notch signaling cascade. Increasing the levels of external Delta ligand (D_ext) can drive EMT through Notch signaling mediated activation of SNAIL, subsequently affecting the dynamics of EMT and PAGE4-AR circuits. We have previously shown that the dynamics of EMT circuit remain largely unaltered upon including miR-34 [47]; thus, we included miR-34 in our circuit, given its connections with Notch signaling circuit [57] (Figure 6A).

The EMT-Notch circuit (i.e., without incorporating the PAGE4-AR circuit) can assume up to four distinct phenotypes when induced with varying levels of external Delta ligand (D_ext) (Figure 6B, Figure S5). In particular, when comparing to the standalone EMT bifurcation (see Figure 1(Bi)), we notice that the epithelial branch (high miR-200) splits into two distinct branches in the EMT-Notch coupled circuit: a (high miR-200, high Notch), i.e., epithelial Receiver state ((E), (R)), and a (high miR-200, high Delta), i.e., epithelial Sender state ((E), (S)). On the other hand, both the hybrid E/M (medium miR-200) and M (low miR-200) branches are characterized by a Sender/Receiver Notch state with (high Notch, high Jagged) ((E/M), (S/R) and (M), (S/R)). For these four phenotypes, the EMT status of a cell is decided by its miR-200 (or ZEB) levels, while the status in terms of S, R, or S/R is defined based on levels of molecules in Notch signaling pathway (Notch, Delta, Jagged).

We first consider the case of strong AR-to-ZEB signalling and weak ZEB-to-AR signalling (${\lambda }_{AtoZ}$ = 0.1, ${\lambda }_{ZtoA}$ = 0.9), thus studying how AR dynamics propagates to EMT and Notch (Figure 6C and Figure S6A). For a low value of D_ext (D_ext = 100), the cell assumes an epithelial phenotype, and Notch receptors and ligands relax to a constant level (Figure 6C, top panel). A higher D_ext level (D_ext = 400) induces a co-existence of two epithelial phenotypes—one of them (epithelial Sender) does not exhibit oscillations in terms of miR-200 and Jagged and the other one (epithelial Receiver) does (Figure 6C, middle panel; Figure S6A). On further increasing Dext (D_ext = 900), the epithelial Sender phenotype is not observed, perhaps because of strong activation of Notch signaling and consequent inhibition of Delta by Notch Intra-Cellular Domain (NICD). In this scenario, the epithelial Receiver phenotype is maintained which continues to exhibit oscillations in both miR-200 and Jagged (Figure 6C, bottom panel; Figure S6A). Overall, a stronger activation of Notch signaling can drive propagation of oscillations seen in PAGE4-AR circuit to the Notch circuit as well through the intermediary EMT circuit (please note that Notch circuit and PAGE4-AR circuits are connected solely through the EMT circuit here).

Next, we consider the opposite case, i.e., weak AR-to-ZEB signaling and strong ZEB-to-AR signaling (${\lambda }_{AtoZ}$ = 0.9, ${\lambda }_{ZtoA}$ = 0.1) and study how Notch multistability affects AR oscillations (Figure 6D and Figure S6B). At low values of Dext (Dext = 100), all initial conditions converge to an epithelial phenotype where ZEB levels are low and thus cannot inhibit the oscillations in AR (Figure 6D, top panel). At intermediate Dext (Dext = 400), AR can either maintain its oscillatory behavior or relax to a low steady state level (Figure 6D, middle panel). On further increasing Dext (Dext = 900), AR loses its oscillations and saturates at a low steady state (Figure 6D, bottom panel). Thus, in case when ZEB inhibits AR strongly, activation of Notch signaling can dampen AR oscillations through increased ZEB levels.

In case of strong inhibition on both sides (${\lambda }_{AtoZ}$ = 0.1, ${\lambda }_{ZtoA}$ = 0.1), we see trends reminiscent of both scenarios discussed above. At low Dext value (Dext = 100), epithelial Sender phenotype is observed, and AR shows oscillations (Figure S7A). As Dext values are increased, Notch signaling is activated that can increase the levels of SNAIL to a value that supports multistability in EMT. Thus, at Dext = 400, we observe 4 states in coupled EMT-Notch-PAGE4 circuit as seen previously in epithelial Sender ((E), (S)), epithelial Receiver ((E), (R)), hybrid E/M Sender/Receiver ((E/M), (S/R)) and Mesenchymal Sender/Receiver ((M), (S/R)) (Figure S7B). Concurrently, AR dynamics shows the co-existence of oscillations and a low steady state value. On increasing Dext further (Dext = 900), the epithelial sender ((E), (S)) and hybrid E/M Sender/Receiver ((E/M), (S/R)) states are not observed, consistent with trends expected of a stronger activation of Notch (Figure S7C).

3. Discussion

Non-genetic heterogeneity is a central theme in cellular decision-making in embryonic development [58]. Its importance in cancer is beginning to be appreciated, with an increasing quantitative and systems level analysis of such heterogeneity to understand the mechanistic underpinnings [59]. Here, we investigate the coupled dynamics of three signaling networks, each of which is capable of generating non-genetic heterogeneity in a cancer cell population. Two of these circuits (EMT, Notch–Delta–Jagged) are multistable [60,61,62,63], and the third one (PAGE4-AR) can exhibit oscillations [45,46]. In our previous work, we had exemplified how the coupled dynamics of both multistable modules (EMT and Notch–Delta–Jagged signaling) at a multi-cell level drove varied spatiotemporal patterns of E, M and hybrid E/M phenotypes [64]. Here, we interrogate the emergent dynamics of coupling of circuits exhibiting multistable and oscillatory behavior and discuss its implications for prostate cancer cells.

Depending on relative strengths of the effect of ZEB1 on AR and vice-versa, we observed that the stand-alone dynamical features of EMT and AR circuits—multistability and oscillations—could percolate to the other circuit, i.e., EMT circuit may show oscillations and/or AR circuits may exhibit bistability. While many features of multistability in EMT has been experimentally observed in multiple cancer cell lines (co-existence of multiple states [65], hysteresis [66,67], and spontaneous state switching [44]), oscillations in EMT have not yet been observed. Encouraging results from preliminary studies using GFP-AR to quantify nuclear translocation of AR in individual PCa cells [68] suggests that empirically validating our theoretical results is technically feasible.

The bidirectional coupling between AR signaling and EMT offer a possible mechanistic link associating the progression of cells towards a partial or full EMT and gain of therapy resistance in prostate cancer. In particular, our model suggests that the epithelial phenotype usually co-occurs with PAGE4 oscillations, whereas transitions to hybrid E/M or mesenchymal phenotypes quench these oscillations and promote low AR levels. From a clinical standpoint, low levels of AR suggest an androgen independent (AI) phenotype that is potentially less susceptible to androgen deprivation therapies. Therefore, EMT induction can potentially promote therapy resistance by stabilizing an androgen independent PCa phenotype through the ZEB1-AR signaling axis. A partial and/or full EMT has been associated with resistance to various therapies in multiple cancers [69]; however, a comparative analysis of partial and full EMT in terms of therapeutic resistance is beyond the scope of the model considered here. Nonetheless, such AR-ZEB1 coupling suggests that not only EMT can drive the acquisition of a drug-resistant state as observed in vitro and in vivo [70,71], but conversely, a switch from drug-sensitive to drug-resistant state can also trigger EMT. Residual breast cancers after conventional therapies have been shown to be mesenchymal [72], but this analysis was at a bulk population level. Thus, these results preclude us from discerning whether the residual cells are a result of phenotypic switching or they are derived from pre-existing mesenchymal cells that were selected during the therapy.

Our model predicts that besides the EMT circuit, Notch signaling pathway can also exhibit oscillations. This pathway has been shown to display oscillatory dynamics, but mostly in the context of somite segmentation clock [50]. The oscillations predicted by our model for EMT and Notch signaling is predicated on a) oscillations in AR levels that remain to be experimentally tested, and b) the mutual inhibition between ZEB1 and AR. In vivo imaging has revealed oscillatory dynamics of ER (estrogen receptor) transcriptional activity in a tissue-dependent manner [73]. Negative feedback loops such as p53-MDM2 can give rise to oscillations [54]; whether such a loop exists for ER and AR (in addition to the one including PAGE4) in PCa cells, remains to be determined. Furthermore, the mutual inhibition between ZEB1 and AR may only be true for PCa cells, but not for other cancers [74]. Any direct coupling between PAGE4/AR and Notch–Delta–Jagged signaling can also alter the emergent dynamics shown for these circuits. Finally, the dynamics of this coupled circuit can be altered by stochasticity in conformational space which can modulate the likelihood of interactions constituting this network. For instance, higher noise in protein conformation can lead to more promiscuous protein interactions, leading to “dynamic networks”, i.e network topology changing as a function of time. Such dynamics can affect the timescale of oscillations and/or the mean residence time in a given stable state (phenotype).

Considered together, our work highlights how coupling between Notch, EMT and PAGE4/AR signaling pathways can give rise to non-trivial dynamics and non-genetic heterogeneity in a cancer cell population. Thus, our results add to the growing literature on investigating coupled dynamics of EMT and other regulatory modules [75,76,77,78], and are likely to impact how we treat PCa in the future.

4. Materials and Methods

Mathematical Modelling of the Coupled PAGE4/AR-EMT-Notch Circuit

In our model, the temporal dynamics of micro-RNAs or proteins (say, X) obey the chemical rate equations that describe their interactions with other species in the circuit. The generic form of such equations is:

$$\frac{dX}{dt}={\Gamma }_X{H}^S\left(A\left(t\right),X\right)-{\gamma }_{X}X$$

(1)

Here, the first term in RHS (${\Gamma }_X{H}^S\left(A\left(t\right),X\right)$), stands for the net production rate of X. This production rate depends on the basal production rate of X (${\Gamma }_X$) and a function $H^S\left(A\left(t\right),X\right)$ which represents the regulation of levels of X resulting from interactions with any another species in the circuit (say, A). In case of regulation on X from multiple species, those terms are multiplied. The second term on the RHS ($-{\gamma }_{X}X$), represents the first-order degradation of X. We used Shifted Hill function to model the effect of one species on another, which leads Equation (1) to take the following form:

$$\frac{dX}{dt}={\Gamma }_X{H}^S\left(A\left(t\right),A0X,n_{AtoX},{\lambda }_{AtoX}\right)-{\gamma }_{X}X$$

(2)

where H^S takes the following form:

$$H^S\left(A\left(t\right),A0X,n_{AtoX},{\lambda }_{AtoX}\right)= \frac{1}{1+{\left(\frac{A\left(t\right)}{A0X}\right)}^{n_{AtoX}}}+{\lambda }_{AtoX}\frac{{\left(\frac{A\left(t\right)}{A0X}\right)}^{n_{AtoX}}}{1+{\left(\frac{A\left(t\right)}{A0X}\right)}^{n_{AtoX}}}$$

The first argument in the Shifted Hill function (A(t)) is number of molecules of the species A at any given time (t). Furthermore, in the PAGE4 circuit, we consider that case of delay in the inhibition of CLK2 by AR. In this particular case, the levels of CLK2 at time t depend on the levels of AR at $\left(t-\tau \right),$ where $\tau$ expresses the amount of delay. A0X is the half-maximal concentration parameter of the Hill function. ${\lambda }_{AtoX}$ is the fold-change of X due to the effect of A. ${\lambda }_{AtoX} >1$ corresponds to activation, and a higher value of ${\lambda }_{AtoX}$ implies a stronger activation. Conversely, ${\lambda }_{AtoX} <1$ corresponds to inhibition, and the closer the value of ${\lambda }_{AtoX}$ to 0, the stronger the inhibition strength. ${\lambda }_{AtoX} =1$ corresponds to the limit case where X does not affect A. Additionally, $n_{AtoX}$ is the characteristic Hill function coefficient which determines the steepness of the Hill function.

Detailed equations for all species in the PAGE4, EMT and Notch circuits are presented in the Supplementary Information (Sections S1–Section S3).