An Effective Model of the Retinoic Acid Induced HL-60 Differentiation Program

In this study, we present an effective model All-Trans Retinoic Acid (ATRA)-induced differentiation of HL-60 cells. The model describes reinforcing feedback between an ATRA-inducible signalsome complex involving many proteins including Vav1, a guanine nucleotide exchange factor, and the activation of the mitogen activated protein kinase (MAPK) cascade. We decomposed the effective model into three modules; a signal initiation module that sensed and transformed an ATRA signal into program activation signals; a signal integration module that controlled the expression of upstream transcription factors; and a phenotype module which encoded the expression of functional differentiation markers from the ATRA-inducible transcription factors. We identified an ensemble of effective model parameters using measurements taken from ATRA-induced HL-60 cells. Using these parameters, model analysis predicted that MAPK activation was bistable as a function of ATRA exposure. Conformational experiments supported ATRA-induced bistability. Additionally, the model captured intermediate and phenotypic gene expression data. Knockout analysis suggested Gfi-1 and PPARg were critical to the ATRAinduced differentiation program. These findings, combined with other literature evidence, suggested that reinforcing feedback is central to hyperactive signaling in a diversity of cell fate programs.


Results
We constructed an effective model of ATRA-induced HL-60 differentiation which described signaling and gene expression events following the addition of ATRA (Fig. 1). HL-60 is a NCI-60 cell line that is a widely used model for studying the mechanism of action of ATRA. There is an immense body of literature on HL-60, therefore it was chosen for modeling. For molecules where there was not explicit HL-60 literature, then literature for well-known interactions from other experimental systems was curated to fill lacuna for the modeling. The model connectivity was developed from literature and the studies presented here (Table 1). We decomposed the ATRA program into three modules; a signal initiation module that sensed and transformed the ATRA signal into activated cRaf-pS621 and the ATRA-RAR/RXR (Trigger) signals (Fig. 1A); a signal integration module that controlled the expression of upstream transcription factors given cRaf-pS621 and activated Trigger signals (Fig. 1B); and a phenotype module which encoded the expression of functional differentiation markers from the ATRA-inducible transcription factors (Fig. 1C). In particular, Trigger (a surrogate for the RARα/RXR transcriptional complex) regulated the expression of the transcription factors CCATT/enhancer binding protein α (C/EBP α), PU.1, and Egr-1. In turn, these transcription factors, in combination with cRaf-pS621, regulated the expression of downstream phenotypic markers such as CD38, CD11b or p47Phox. Each component of these modules was described by a mRNA and protein balance equation. Additionally, the signal initiation module also described the abundance of activated species e.g., Trigger and cRaf-pS621 whose values were derived from unactivated Trigger and cRaf protein levels. Lastly, because the population of HL-60 cells was dividing, we also considered a dilution term in all balance equations. The signal initiation module contained nine differential equations, while the signal integration and phenotype modules were collectively encoded by 54 differential equations. Model parameters were taken from literature (Table 2), or estimated from experimental data using heuristic optimization (see materials and methods).
The signal initiation module recapitulated sustained signalsome and MAPK activation following exposure to 1 μM ATRA ( Fig. 2A,B). An ensemble of effective model parameters was estimated by minimizing the difference between simulations and time-series measurements of BLR1 mRNA and cRaf-pS621 following the addition of 1 μM ATRA. We focused on the S621 phosphorylation site of cRaf since enhanced phosphorylation at this site is a defining characteristic of sustained MAPK signaling activation in HL-60. The effective model captured both ATRA-induced BLR1 expression ( Fig. 2A) and sustained phosphorylation of cRaf-pS621 (Fig. 2B) in a growing population of HL-60 cells. Together, the reinforcing feedback within the signalsome and its embedded MAPK signaling axis led to sustained activation over multiple cellular generations. However, the effective model failed to capture the decline of BLR1 message after 48 hr of ATRA exposure. This suggested that we captured the logic leading to the onset of differentiation, but failed to describe program shutdown. Much of the focus in the literature has been on understanding the initiation of differentiation, with little attention paid to understanding how a program is terminated. This is a potential new direction that could be explored. Next, we tested the response of the signal initiation module to different ATRA dosages. The signal initiation model was bistable with respect to ATRA induction (Fig. 2C,D). Phaseplane analysis predicted two stable steady-states when ATRA was present below a critical threshold (Fig. 2C), and only a single steady-state above the threshold (Fig. 2D). In the lower stable state, neither the signalsome nor cRaf-pS621 were present (thus, the differentiation program was inactive). However, at the higher stable state, both the signalsome and cRaf-pS621 were present, allowing for sustained activation and differentiation. Interestingly, when ATRA was above a critical threshold, only the activated state was accessible (Fig. 2D). To test these findings, we first identified the ATRA threshold. We exposed HL-60 cells to different ATRA concentrations for 72 hr (Fig. 2E). Morphological changes associated with differentiation were visible for ATRA ≥0.25 μM, suggesting the critical ATRA threshold was near this concentration. Next, we conducted ATRA washout experiments to determine if activated cells remained activated in the absence of ATRA. HL-60 cells locked into an activated state remained activated following ATRA withdraw (Fig. 3C). This sustained activation resulted from reinforcing feedback between the signalsome and the MAPK pathway. Thus, following activation, if we inhibited or removed elements from the signal initiation module we expected the signalsome and MAPK signals to decay. We simulated ATRA induced activation in the presence of kinase inhibitors, and without key circuit elements. Consistent with experimental results using multiple MAPK inhibitors, ATRA activation in the presence of MAPK inhibitors lowered the steady-state value of signalsome (Fig. 3A). In the presence of BLR1, the signalsome and cRaf-pS621 signals were maintained following ATRA withdraw (Fig. 3B, gray). On the other hand, BLR1 deletion removed the ability of the circuit to maintain a sustained MAPK response following the withdraw of ATRA (Fig. 3B, blue). Lastly, washout experiments in which cells were exposed to 1 μM ATRA for 24 hr, and then transferred to fresh media without ATRA, confirmed the persistence of the self sustaining activated state for up to 144 hr (Fig. 3C). Thus, these experiments confirmed that reinforcing positive feedback likely drives the ATRA-induced differentiation program. Next, we analyzed the ATRA-induced downstream gene expression program following signalsome and cRaf activation.  The signal integration and phenotype modules described ATRA-induced gene expression in wild-type HL-60 cells (Fig. 4). The signal initiation module produced two outputs, activated Trigger and cRaf-pS621 which drove the expression of ATRA-induced transcription factors, which then in turn activated the phenotypic program. We assembled the connectivity of the signal integration and phenotypic programs driven by Trigger and cRaf-pS621 from literature (Table 1). We estimated the parameters for the signal initiation, and phenotype modules from steady-state and dynamic measurements of transcription factor and phenotypic marker expression following the addition of ATRA [28][29][30][31] . However, the bulk of the model parameters were taken from literature 32 and were not estimated in this study (see materials and methods). The model simulations captured the time dependent expression of CD38 and CD11b following the addition ATRA (Fig. 4A), and the steady-state for signal integration and phenotypic markers (Fig. 4B). Lastly, we used the predicted values of the p21 and E2F protein abundance to estimate a blackbox model of ATRA-induced G0 arrest (Fig. 5). The phenotype module predicted p21 expression significantly increased and E2F expression decreased, in response to ATRA exposure (Fig. 5A). We then used the ratio of these values in a polynomial model to calculate the fraction of HL-60 cells in G0 arrest following the addition of ATRA (Fig. 5B). The third-order polynomial model captured the trend in measured G0-arrest values as a function of time, and was robust to uncertainty in the measured data (Fig. 5B, gray). Taken together, the output of the signal integration and phenotypic modules was consistent with time-series and steady-state measurements, thereby validating the assumed molecular connectivity. Moreover, outputs from the phenotype module described the trend in ATRA-induced G0 cell cycle arrest. Next, we explored which proteins and protein interactions in the signal integration module most influenced the system response.
The Gfi-1 and PPARγ proteins were important regulators of ATRA-induced signal integration and phenotypic change (Fig. 6). We conducted pairwise gene knockout simulations in the signal integration and phenotype modules to estimate which proteins controlled the processing of the Trigger and cRaf-S621 signals. The difference between the system state with and without the gene knockouts (encoded as a normalized state displacement matrix) was decomposed using Singular Value Decomposition (SVD). A panel of ten parameter sets was sampled, and the average normalized displacement matrix was decomposed. The first six modes (approximately 36% of the total) described ≥95% of the gene knockout variance, with the most important components of these modes being the Gfi-1 and PPARγ proteins, and to a lesser extent PU.1, C/EBPα and and AP1 (Fig. 6A). To better understand which protein-DNA connections were important, we simulated the pairwise deletion of interactions between these proteins and their respective regulatory targets. Singular value decomposition of the normalized state displacement matrix assembled from the pairwise connection deletions, suggested the first six modes (approximately 26% of the total) accounted for ≥90% of the variance. Globally, the most sensitive interactions controlled p47Phox and p21 expression, markers for the cell-cycle arrest and reactive oxygen phenotypic axes activated following ATRA addition (Fig. 6B). While the p21 spot appeared small, it was the second highest  ranked response behind p47Phox, in the largest response mode. The interactions associated with these shifts likely involved important components; the deleted interactions involved the action of PU.1, C/EBPα and cRaf at both the p47Phox and p21 promoters, as well as PPARγ action for p21. Taken together, the gene and interaction knockout studies showed that the action of PPARγ, Gfi-1 and C/EBPα was consistently important over multiple target genes. The connection knockout analysis also revealed robustness within the network. For example, no pair of deletions qualitatively changed the expression of regulators such as PU.1, Oct1, Oct4 or PPARγ. Thus, the expression of these species was robust to disturbance in the connectivity. To better understand the combined influence of the PPARγ and Gfi-1 deletions, we computed the fold change in the protein levels in the single (Gfi-1 −/− or PPARγ −/− ) and double (Gfi-1 −/− and PPARγ −/− ) mutants for the best fit parameter set (Fig. 7). Deletion of Gfi-1 led to a 2-4 fold increase in EGR-1, CD11b and C/EBPα expression, and a >8 fold increase in PU.1 abundance ( Fig. 7,blue). On the other hand, deletion of PPARγ led to >8 fold downregulation of CD38, p21, IRF1 and Oct1 (Fig. 7, red). Both knockouts slightly increased E2F expression, but neither influenced the expression of p47Phox. The double mutant was qualitatively similar to the combined behavior of the two single mutant cases. Taken together, Gfi-1 and PPARγ controlled the cell-cycle arrest and receptor signaling axes, with PPARγ regulating CD38, IRF1 and p21 expression while Gfi-1 controlled CD11b expression. These simulations suggested deletion of PPARγ and Gfi-1 would not interfere with reactive oxygen formation, but would limit the ability of HL-60 cells to arrest. However, this analysis did not give insight into which components upstream of the signal initiation module were important. Toward this question, we explored the composition and regulation of the signalsome complex by experimentally interrogating a panel of possible Raf interaction partners. The full composition of the signalsome, and the kinase therein ultimately responsible for mediating ATRA-induced Raf activation is still not known. To explore this question, we conducted immunoprecipitation and subsequent Western blotting to identify interactions between Raf and 19 putative interaction partners. A panel of 19 possible Raf interaction partners (kinases, GTPases, scaffolding proteins etc) was constructed based upon known signaling pathways. We did not consider the most likely binding partner, the small GTPase RAS, as previous studies have ruled it out for MAPK activation in HL-60 cells 25,33 . Total Raf was used as a bait protein for the immunoprecipitation studies. Interrogation of the Raf interactome thus suggested Vav1 was involved with ATRA-induced initiation of MAPK activity (Fig. 8). Western blot analysis using total Raf and Raf-pS621 specific antibodies confirmed the presence of the bait protein, total and phosphorylated forms, in the immunoprecipitate (Fig. 8A). Of the 19 proteins sampled, Vav1, Src, CK2, Akt, and 14-3-3 co-precipitated with Raf, suggesting their co-existence in a complex was possible. However, only the associations between Raf and Vav1, and Raf and Src were ATRA-inducible (Fig. 8). The interaction between Vav1 and Raf was one of the most prominent interactions in the panel, and it was crippled by inhibiting Raf. Furthermore, the Vav1 and Src associations were correlated with Raf-pS621 abundance in the precipitate. Other proteins e.g., CK2, Akt and 14-3-3, generally bound Raf regardless of phosphorylation status or ATRA treatment. The remaining proteins sampled were expressed in whole cell lysate (Fig. 8B), but were not detectable in the immuno-precipitate with Raf IP; consistent with the potential importance of the Raf-Vav interaction for signaling, it paralleled Raf phosphorylation at S621, a putative telltale of the activated kinase. Furthermore, treatment with the Raf kinase inhibitor GW5074 following ATRA exposure reduced the association of both Vav1 with Raf and Src with Raf ( Fig. 8C), although the signal intensity for Src was notably weak. However, GW5074 did not influence the association of CK2 or 14-3-3 with Raf, further demonstrating their independence from Raf phosphorylation. Interestingly, the Raf-Akt interaction qualitatively increased following treatment with GW5074; however, it remained unaffected by treatment with ATRA. Src family kinases are known to be important in myeloid differentiation 34 and their role in HL-60 differentiation has been  Gene expression at t = 48 hr following ATRA exposure. Gene expression was normalized to expression in the absence of ATRA. The gene expression is quantified by the protein fold change of quantified Western blot data (from at least three biological repeat nuclear lysates) using ImageJ. Experimental data in panels A and B were reproduced from Jensen et al. 31 . Model simulations were conducted using the ten best parameter sets collected during model identification. Solid lines (or bars) denote the mean model performance, while the shaded region (or error bars) denote the 95% confidence interval calculated over the parameter ensemble. investigated elsewhere 15 . Given the existing work and variable reproducibility in the context of the Raf immunoprecipitate, we did not investigate the role of Src further in this study. Taken together, the immunoprecipitation and GW5074 results implicated Vav1 association to be correlated with Raf activation following ATRA-treatment. Singular value decomposition of the average system response (l 2 -norm between the perturbed and nominal state) following pairwise gene knockout simulations using the top ten best fit parameter sets. The rows denote the deleted genes, while columns denote the response mode. (B) Singular value decomposition of the average system response (l 2 -norm between the perturbed and nominal state) following the pairwise removal of protein-DNA connections for the top ten best fit parameter sets. The rows denote protein-DNA interactions at the labeled promoter, while the columns denote the top ranked response modes (combinations of deletions). The percentage at the top of each column describes the fraction of the variance in the system state captured by the node combinations in the rows. Robustness of the HL-60 differentiation program following exposure to 1 μM ATRA at t = 0 hr. Protein fold change at t = 48 hr (rows) in single and double knock-out mutants (columns) relative to wild-type HL-60 cells. The responses were grouped into >2,4 and 8 fold changes. The best fit parameter set was used to calculate the protein fold change.
Further, while we observed possible immunoprecipitation of Src with Raf, the western blot results showed inconsistent results and significant non-specific binding; therefore we could not rule in or out a Src/Raf interaction. Previous studies demonstrated that a Vav1-Slp76-Cbl-CD38 complex plays an important role in ATRA-induced MAPK activation and differentiation of HL-60 cells 17 . Here we did not observe direct interaction of Raf with Cbl or Slp76; however, this interaction could could be involved upstream. Next, we considered the effect of the Raf kinase inhibitor GW5074 on functional markers of ATRA-induced growth arrest and differentiation.
Inhibition of Raf kinase activity modulated MAPK activation and differentiation markers following ATRA exposure (Fig. 8D-F). ATRA treatment alone statistically significantly increased the G1/G0 percentage over the untreated control, while GW5074 alone had a negligible effect on the cell cycle distribution (Fig. 8D). Surprisingly, the combination of GW5074 and ATRA statistically significantly increased the G1/G0 population (82 ± 1%) compared with ATRA alone (61 ± 0.5%). Increased G1/G0 arrest following the combined treatment with GW5074 and ATRA was unexpected, as the combination of ATRA and the MEK inhibitor (PD98059) has been shown previously to decrease ATRA-induced growth arrest 12 . However, growth arrest is not the sole indication of functional differentiation. Expression of the cell surface marker CD11b has also been shown to coincide with HL-60 cells myeloid differentiation 35 . We measured CD11b expression, for the various treatment groups, using immuno-fluorescence flow cytometry 48 hr post-treatment. As with G1/G0 arrest, ATRA alone increased CD11b expression over the untreated control, while GW5074 further enhanced ATRA-induced CD11b expression (Fig. 8E). GW5074 alone had no statistically significant effect on CD11b expression, compared with the untreated control. Lastly, the inducible reactive oxygen species (ROS) response was used as a functional marker of differentiated neutrophils 20 . We measured the ROS response induced by the phorbol ester 12-O-tetradecanoylphorbol-13-acetate (TPA) using flow cytometry. Untreated cells showed no discernible TPA response, with only 7.0 ± 3.0% ROS induction (Fig. 8F). Cells treated with ATRA had a significantly increased TPA response, 53 ± 7% ROS induction 48 hr post-treatment. Treatment with both ATRA and GW5074 statistically significantly reduced ROS induction (22 ± 0.6%) compared to ATRA alone. Interestingly, Western blot analysis did not detect a GW5074 effect on ATRA-induced expression of p47Phox, a required upstream component of the ROS response (Fig. 8F,  bottom). Thus, the inhibitory effect of GW5074 on inducible ROS might occur downstream of p47Phox expression. However, the ROS producing complex is MAPK dependent, therefore it is also possible that GW5074 inhibited ROS production by interfering with MAPK activation (in which case the p47Phox marker might not accurately reflect phenotypic conversion and differentiation).

Discussion
In this study, we presented an effective model of ATRA-inducible differentiation of HL-60 cells. The model consisted of three modules: a signal initiation module that sensed and transformed the ATRA signal into activated cRaf-pS621 and the ATRA-RAR/RXR (Trigger) signals; a signal integration module that controlled the expression of upstream transcription factors given cRaf-pS621 and activated Trigger signals; and a phenotype module which encoded the expression of functional differentiation markers from the ATRA-inducible transcription factors. The model described the transcription and translation of genes in each module, and signaling events in each module in a growing population of HL-60 cells. Model parameters were taken from literature, however, unknown coefficients that appear in the promoter logic models were estimated from protein measurements in HL-60 cells following ATRA exposure. Despite its simplicity, the effective model captured key features of the ATRA induced differentiation such as sustained MAPK activation, and bistability with respect to ATRA exposure. The model also described the expression of upstream transcription factors which regulated the expression of differentiation markers. Lastly, analysis of the response of the model to perturbations identified Gfi-1 and PPARγ as master regulators of ATRA-induced differentiation. We also found evidence of a prominent role for an ATRA-inducible component of the signalsome, Vav1. Vav1 is a guanine nucleotide exchange factor for Rho family GTPases that activate pathways leading to actin cytoskeletal rearrangements and transcriptional alterations 36 . The Vav1/Raf association correlated with Raf activity, was ATRA-inducible and decreased after treatment with the Raf inhibitor GW5074.
Naturally occurring cell fate decisions often incorporate reinforcing feedback and bistability 37,38 . One of the most well studied cell fate circuits is the Mos mitogen-activated protein kinase cascade in Xenopus oocytes. This cascade is activated when oocytes are induced by the steroid hormone progesterone 39 . The MEK-dependent activation of p42 MAPK stimulates the accumulation of the Mos oncoprotein, which in turn activates MEK, thereby closing the feedback loop. This is similar to the signal initiation module presented here; ATRA drives signalsome formation, which activates MAPK, which in turn leads to more signalsome activation. Thus, while HL-60 and Xenopus oocytes are vastly different biological models, their cell fate programs share a similar architectural feature. Reinforcing feedback and bistability has also been implicated in hematopoietic cell fate determination. Laslo et al. showed in nonmalignant myelomonocytic cells that the counter antagonistic repressors, Gfi-1 and Egr-1/2 (whose expression is tuned by PU.1 and C/EBPα), encode a bistable switch that results in a macrophage, neutrophil or a mixed lineage population depending upon PU.1 and C/EBPα expression 38 . The current model contained the Gfi-1 and Egr-1/2 agonistic switch; however, its significance was unclear for HL-60 cells. The expression of Gfi-1, Egr-1/2, C/EBPα and PU.1 was not consistent with the canonical lineage pattern expected from literature. For example, Egr-1/2 expression (associated with a macrophage lineage) increased, while Gfi-1 expression (associated with a neutrophil lineage) was unchanged following ATRA exposure. Thus, HL-60 cells, which are a less mature cancer cell line, exhibited a non-canonical expression pattern. Other unrelated cell fate decisions such as programmed cell death have also been suggested to be bistable 40 . Still more biochemical networks important to human health, for example the human coagulation or complement cascades, also feature strong positive feedback elements 41 . Thus, while reinforcing feedback is often undesirable in human engineered systems, it is at the core of a diverse variety of cell fate programs and other networks important to human health.
Analysis of the signal integration and phenotype modules suggested Gfi-1 and PPARγ proteins were important regulators of ATRA-induced signal integration and phenotypic change. Model analysis showed that PU.1, Egr-1 and C/EBPα expression increased in Gfi-1 −/− mutants, where PU.1 expression was upregulated by greater than 8-fold. Simulations suggested that combined Gfi-1 + PPARg deletion crippled the ability of HL-60 cells to undergo neutrophilic differentiation following ATRA exposure. This confirms previous literature showing that Gfi-1 KO mice lack normal neutrophils 42 . PU.1, a member of the ets transcription factor family, is a well known regulator of granulocyte and monocyte development 43 . The relative level of PU.1 and C/EBPα is thought to control macrophage versus neutrophil cell fate decisions in granulocytic macrophage progenitor cells 44 . Simulations suggested that combined Gfi-1 + PPARγ deletion crippled the ability of HL-60 cells to undergo neutrophilic differentiation following ATRA exposure. p21 expression decreased significantly, suggesting Gfi-1 −/− + PPARγ −/− mutants were less likely to G0-arrest following ATRA exposure. The expression of other neutrophilic markers, such as CD38, also decreased in Gfi-1 −/− + PPARγ −/− cells. On the other hand, the expression of reactive oxygen metabolic markers, or other important transcription factors such as Oct4 did not change. For example, model analysis suggested that the C/EBPα dependent interaction of PU.1 with the NCF1 gene, which encodes the p47Phox protein, was the most sensitive PU.1 connection; deletion of this connection removed the ability of the system to express p47Phox. p47Phox, also known as neutrophil cytosol factor 1, is one of four cytosolic subunits of the multi-protein NADPH oxidase complex found in neutrophils 45 . This enzyme is responsible for reactive oxygen species (ROS) production, a key component of the anti-microbial function of neutrophils. While p47Phox expression required C/EBPα and PU.1, neither Gfi-1 nor PPARγ deletion increased expression. This suggested that p47Phox expression was saturated with respect to C/EBPα and PU.1, and simultaneously not sensitive to PPARγ abundance. Taken together, Gfi-1 −/− + PPARγ −/− cells were predicted to exhibit some aspects of the ATRA response, but not other critical features such as cell cycle arrest. Hock et al. showed that Gfi-1 −/− mice lacked normal neutrophils, and were highly sensitive to bacterial infection 42 . Thus, the model analysis was consistent with this study. However, other predictions concerning the behavior of the Gfi-1 −/− + PPARγ −/− mutants remain to be tested.
Immunoprecipitation studies identified a limited number of ATRA-dependent and -independent Raf interaction partners. We established potential interactions between Raf and key partners such as Vav1, Src, Akt, CK2 and 14-3-3. However, we were unable to detect the association of Raf with common kinases and GTPases such as PKC, PKA, p38, Rac and Rho as observed in literature 46,47 . To investigate the association of c-Raf and PKA or PKC, we assessed both the expression levels and the associations of those important signaling molecules. Surprisingly, the expression levels of PKA, PKCα and PKCγ were not ATRA regulated, nor did we detect an association with Raf. We believe that our current perception of the signaling pathway driving differentiation in this model is novel. It also diverges from the classical perception in that its activation is not Ras driven. This has been reported by Katagiri et al. that Ras is not a driver in RA-induced differentiation in HL-60 33 . Finally, we note that the classical paradigms were typically derived in NIH3T3 cells where signal duration is a rapid transient. By contrast in HL-60, it is a prolonged MAPK signaling that drives RA-induced granulocytic differentiation.
All of these partners are known to be associated with Raf activation or function. Src is known to bind Raf through an SH2 domain, and this association has been shown to be dependent of the serine phosphorylation of Raf 48 . Thus, an ATRA inducible Src/Raf association may be a result of ATRA-induced Raf phosphorylation at S259 or S621. We also identified an interaction between Raf and the Ser/Thr kinases Akt and CK2. Akt can phosphorylate Raf at S259, as demonstrated by studies in a human breast cancer line 49 . CK2 can also phosphorylate Raf, although the literature has traditionally focused on S338 and not S621 or S259 50 . However, neither of these kinase interactions were ATRA-inducible, suggesting their association with Raf alone was not associated with ATRA-induced Raf phosphorylation. The adapter protein 14-3-3 was also constitutively associated with Raf. The interaction between Raf and 14-3-3 has been associated with both S621 and S259 phosphorylation and activity 51 . Additionally, the association of Raf with 14-3-3 not only stabilized S621 phosphorylation, but also reversed the S621 phosphorylation from inhibitory to activating 52 . Finally, we found that Vav1/Raf association correlated with Raf activity, was ATRA-inducible and decreased after treatment with GW5074. The presence of Vav1 in Raf/ Grb2 complexes has been shown to correlate with increased Raf activity in mast cells 53 . Furthermore, studies on Vav1 knockout mice demonstrated that the loss of Vav1 resulted in deficiencies of ERK signaling for both T-cells as well as neutrophils 54,55 . Interestingly, while an integrin ligand-induced ROS response was blocked in Vav1 knockout neutrophils, TPA was able to bypass the Vav1 requirement and stimulate both ERK phosphorylation and ROS induction 55 . It is possible that Vav1 is downstream of various integrin receptors but upstream of Raf in terms of inducible ROS responses. Vav1 has also been shown to associate with a Cbl-Slp76-CD38 complex in an ATRA-dependent manner; furthermore, transfection of HL-60 cells with Cbl mutants that fail to bind CD38, yet still bind Slp76 and Vav1, prevents ATRA-induced MAPK activation 17 . The literature suggest a variety of possible receptor-signaling pathways, which involve Vav1, for MAPK activation; moreover, given the ATRA-inducible association Vav1 may play a direct role in Raf activation.
We hypothesized that Vav1 is a member of an ATRA-inducible signalsome complex which propels sustained MAPK activation, arrest and differentiation (shown schematically in Fig. 9). Initially, ATRA-induced Vav1 expression drives increased association between Vav1 and Raf. This increased interaction facilitates phosphorylation and activation of Raf by pre-bound Akt and/or CK2 at S621 or perhaps S259. Constitutively bound 14-3-3 may also stabilize the S621 phosphorylation, modulate the activity and/or up-regulate autophosphorylation. Activated Raf can then drive ERK activation, which in turn closes the positive feedback loop by activating Raf transcription factors e.g., Sp1 and/or STAT1 [56][57][58][59] . We tested this working hypothesis using mathematical modeling. The model recapitulated both ATRA time-course data as well as the GW5074 inhibitor effects. This suggested the proposed Raf-Vav1 architecture was at least consistent with the experimental studies. Further, analysis of the Raf-Vav1 model identified bistability in phosphorylated ERK levels. Thus, two possible MAPK activation branches were possible for experimentally testable ATRA values. The analysis also suggested the ATRA-induced Raf-Vav1 architecture could be locked into a sustained signaling mode (high phosphorylated ERK) even in the absence of a ATRA signal. This locked-in property could give rise to an ATRA-induction memory. We validated the treatment memory property predicted by the Raf-Vav1 circuit experimentally using ATRA-washout experiments. ERK phosphorylation levels remained high for more then 96 hr after ATRA was removed. Previous studies demonstrated that HL-60 cells possessed an inheritable memory of ATRA stimulus 60 . Although the active state was self-sustaining, the inactive state demonstrated considerable robustness to perturbation. For example, we found that 50x overexpression of Raf was required to reliably lock MAPK into the activated state, while small perturbations had almost no effect on phosphorylated ERK levels over the entire ensemble. CD38 expression correlated with the phosphorylated ERK, suggesting its involvement in the signaling complex. Our computational and experimental results showed that positive feedback, through ERK-dependent Raf expression, could sustain MAPK signaling through many division cycles. Such molecular mechanisms could underly aspects of cellular memory associated to consecutive ATRA treatments.

Methods
Effective gene expression model equations. The ATRA differentiation model was encoded as a system of differential algebraic equations (DAEs) which described both signaling and gene expression processes. We modeled transcription and translation as Ordinary Differential Equations (ODEs), while signaling processes were assumed to quickly equilibrate and were treated as a pseudo steady state system of algebraic equations. The model formulation follows from a previous study of the Epithelial Messecnchymal Transition (EMT) 61 ; in the current study additional attention was paid to the formulation of the transcription and translation rates, and an updated approach was taken to model the regulation of gene expression.
We decomposed the ATRA-induced differentiation program into three modules; a signal initiation module that sensed and transformed the ATRA signal into activated cRaf-pS621 and the ATRA-RAR/RXR (activated Trigger) signals; a signal integration module that controlled the expression of upstream transcription factors given cRaf-pS621 and activated Trigger signals; and a phenotype module which encoded the expression of functional differentiation markers from the ATRA-inducible transcription factors. The output of the signal initiation module was the input to the gene expression model. For each gene j 1,2, ,  = … , we modeled both the mRNA (m j ), protein (p j ) and signaling species abundance: Figure 9. This schematic diagram shows the hypothetical principal pathways in the ATRA-induced signaling that results in cell differentiation in the HL-60 myeloid leukemia model 17,[67][68][69][70][71] . It is based on modules and feedback loops. There are three main arms (shown top to bottom): 1. Direct ATRA targeting of RAREs in genes such as CD38 or BLR1; 2. Formation of a signalsome that has a regulatory module that includes Vav (a guanine nucleotide exchange factor), CBL and SLP-76 (adaptors), and Lyn (a Src family kinase) that regulates a Raf/Mek/ Erk axis that incorporates Erk to Raf feedback, where the regulators are modulated by AhR and CD38 receptors; and 3. Direct ATRA targeted up regulation of CDKI to control RB hypophosphorylation. The Raf/Mek/Erk axis is embedded in the signalsome and subject to modulation by the regulators. The output of the signalsome is discharge of the Raf from the cytosol to the nucleus where it binds (hyper)phospho-RB and other targets, including NFATc3, which enables activation of the ATRA bound RAR/RXR poised on the BLR1 promoter, and also GSK3, phosphorylation of which relieves its inhibitory effect on RAR α. CDKI directed hypophosphorylation of RB releases Raf sequestered by RB to go to NFATc3, GSK3, and other targets. A significant consequence of the nuclear RAF is ergo ultimately to enable or hyperactivate transcriptional activation by RAR α to drive differentiation. It might be noted that this proposed general model provides a mechanistic rationalization for why cell cycle arrest is historically oft times perceived as a precondition for phenotypic maturation.
SCIENtIFIC REPORTS | 7: 14327 | DOI:10.1038/s41598-017-14523-5 . The model parameter vector is denoted by κ. The terms r T j , and r X j , denote the specific rates of transcription, and translation while the terms θ m j , and p j , θ denote first-order degradation constants for mRNA and protein, respectively. The specific transcription rate r T j , was modeled as the product of a kinetic term r T j , and a control term u j which described how the abundance of transcription factors, or other regulators influenced the expression of gene j.
The gene expression control term u 0 1 j ≤ ≤ depended upon the combination of factors which influenced the expression of gene j. If the expression of gene j was influenced by … m 1, , factors, we modeled this relationship as ij denotes a regulatory transfer function quantifying the influence of factor i on the expression of gene j, and  ( ) j ⋅ denotes an integration rule which combines the individual regulatory inputs for gene j into a single control term. In this study, the integration rule governing gene expression was the weighted fraction of promoter configurations that resulted in gene expression 62 : The numerator, the weighted sum (with weights W nj ) of promoter configurations leading to gene expression, was normalized by all possible promoter configurations (denominator). The likelihood of each configuration was quantified by the transfer function f nj (which we modeled using Hill functions), while the lead term in the numerator W R j , 1 denotes the weight of constitutive expression for gene j. Given the formulation of the control law, the j λ term (which denotes the constitutive rate of expression of gene j) was given by: The kinetic transcription term r T j , was modeled as: where the maximum gene expression rate V T max was defined as the product of a characteristic transcription rate constant (k T ) and the abundance of RNA polymerase (R 1 ), , term denotes the ratio of transcription read lengths; L T o , represents a characteristic gene length, while L T j , denotes the length of gene j.
Thus, the ratio ( ) , is a gene specific correction to the characteristic transcription rate V T max . If a gene expression process had no modifying factors, u 1 j = . Lastly, the specific translation rate was modeled as: where V X max denotes a characteristic maximum translation rate estimated from literature, and K X denotes a translation saturation constant. The characteristic maximum translation rate was defined as the product of a characteristic translation rate constant (k X ) and the Ribosome abundance (R 2 ), V k R ( ) As was the case for transcription, we corrected the characteristic translation rate by the ratio of the length of a characteristic transcript normalized by the length of transcript j. The sequence lengths used in this study are given in Table 3; the characteristic gene and mRNA lengths were given by the average lengths computed from the values in Table 3.
Signaling model equations. The signal initiation and integration modules required the abundance of cRaf-pS621 and ATRA-RAR/RXR (activated Trigger) as inputs. However, the base model described only the abundance of inactive proteins e.g., cRaf or RAR/RXR but not the activated forms. To address this issue, we estimated pseudo steady state approximations for the abundance of cRaf-pS621 and activated Trigger. The abundance of activated trigger (x a,1 ) was estimated directly from the RAR/RXR abundance (x u,1 ): where α denotes a gain parameter; α = .
0 0 if ATRA is less than a threshold, and 0 1 α = . if ATRA is greater than the differentiation threshold. The abundance of cRaf-pS621 was estimated by making the pseudo steady state approximation on the cRaf-pS621 balance. In general, the abundance of an activated signaling species i was governed by: The quantity x i denotes concentration of signaling species i, while  and  denote the number of signaling reactions and signaling species in the model, respectively. The term + r x k ( , ) i , denotes the rate of generation of activated species i, while μ denotes the specific growth rate, and k d i , denotes the rate constant controlling the non-specific degradation of x i . We neglected deactivation reactions e.g., phosphatase activities. We assumed that signaling processes were fast compared to gene expression; this allowed us to approximate the signaling balance as: The generation rate was written as the product of a kinetic term ( + r i , ) and a control term (v i ). The control terms v 0 1 j ≤ ≤ depended upon the combination of factors which influenced rate process j. If rate j was influenced by m 1, , … factors, we modeled this relationship as  Estimation of gene expression model parameters. Parameters appearing in the mRNA and protein balances, e.g., maximum transcription and translation rates, the half-life of a typical mRNA and proteins (assumed to be same for all transcripts/proteins), and typical values for the copies per cell of RNA polymerase and ribosomes were estimated from literature ( Table 2). The saturation constants K X and K T appearing in the transcription and translation rate equations were adjusted so that gene expression and translation resulted in gene products on a biologically realistic concentration scale. Lastly, we calculated the concentration for gene G j by assuming, on average, that a cell had two copies of each gene at any given time. Thus, the bulk of our model parameters were taken from literature, and were not adjusted during model identification. However, the remaining parameters, e.g., the W ij values or parameters appearing in the transfer functions f dj which appeared in the gene expression control laws, were estimated from the experimental data discussed here. We assumed promoter configuration weights were bounded between ∈ W [0, 100] ij ; all cooperativity coefficients η ij appearing in the binding transfer functions f dj were bounded between η ∈ [0,4] ij ; and all disassociation constants K ij appearing in the binding transfer functions f dj were bounded between K [0,1000] ij ∈ (nM). Signaling and gene expression model parameters were estimated by minimizing the squared difference between simulations and experimental protein data set j. We measured the squared difference in the scale, fold change and shape for protein j: The first term in Eq. (13) quantified the initial scale error, directly before the addition of ATRA. In this case,  − t ( ) j (the approximate concentration of protein j before the addition of ATRA) was estimated from literature. This term was required because the protein measurements were reported as the fold-change; thus, the data was normalized by a control value measured before the addition of ATRA. However, the model operated on a physical scale. The first term allowed the model to capture physically realistic changes following ATRA addition. The second term quantified the difference in the fold-change of protein j as a function of time. The terms ij  and y iĵ denote the scaled experimental observations and simulation outputs (fold-change; protein normalized by control value directly before ATRA addition) at time i from protein j, where  j denoted the number of time points for data set j. Lastly, the third term of the objective function measured the difference in the shape of the measured and simulated protein levels. The scaled value  0 1 ij ≤ ≤ ′ was given by: and  = ′ 1 ij describe the lowest (highest) intensity bands. A similar scaling was used for the simulation output. We minimized the total model residual E j j ∑ using a heuristic direct-search optimization procedure, subject to box constraints on the parameter values, starting from a random initial parameter guess. Each downhill step was archived and used for ensemble calculations. The optimization procedure (a covariance matrix adaptation evolution strategy) has been reported previously 64 . Estimation of an effective cell cycle arrest model. We formulated an effective N-order polynomial model of the fraction of cells undergoing ATRA-induced cell cycle arrest at time t, t ( )  , as: denotes a basis function. The basis functions were dependent upon the system state; in this study, we used N = 4 and basis functions of the form: The parameters a a , , 0 3 … were estimated directly from cell-cycle measurements (biological replicates) using least-squares. The form of the basis function assumed p21 was directly proportional, and E2F inversely proportional, to G0-arrest. However, this was one of many possible forms for the basis functions.
Cell culture and treatment. Human myeloblastic leukemia cells (HL-60 cells) were grown in a humidified atmosphere of 5% CO 2 at 37 °C and maintained in RPMI 1640 from Gibco (Carlsbad, CA) supplemented with 5% heat inactivated fetal bovine serum from Hyclone (Logan, UT) and 1× antibiotic/antimicotic (Gibco, Carlsbad, CA). Cells were cultured in constant exponential growth 65 . Experimental cultures were initiated at 0 1 10 6 . × cells/mL 24 hr prior to ATRA treatment; if indicated, cells were also treated with GW5074 (2 μM) 18 hr before ATRA treatment. For the cell culture washout experiments, cells were treated with ATRA for 24 hr, washed 3x with prewarmed serum supplemented culture medium to remove ATRA, and reseeded in ATRA-free media as described. Western blot analysis was performed at incremental time points after removal of ATRA.
Chemicals. All-Trans Retinoic Acid (ATRA) from Sigma-Aldrich (St. Louis, MO) was dissolved in 100% ethanol with a stock concentration of 5 mM, and used at a final concentration of 1 μM (unless otherwise noted). The