Introduction

Differentiation induction chemotherapy (DIC), using agents such as the vitamin A derivative all-trans retinoic acid (ATRA), is a promising approach for the treatment of many cancers1,2,3,4,5. For example, ATRA treatment induces remission in 80–90% of promyelocytic leukemia (APL) PML-RARα-positive patients6, thereby transforming a fatal diagnosis into a manageable disease7. However, remission is sometimes not durable and relapsed cases exhibit emergent ATRA resistance8,9. To understand the basis of this resistance, we must first understand the ATRA-induced differentiation program. Toward this challenge, lessons learned in model systems, such as the lineage-uncommitted human myeloblastic cell line HL-60 reported to closely resemble patient derived cells10, could inform our analysis of the differentiation programs occurring in patients. Patient derived HL-60 leukemia cells have been a durable experimental model since the 1970’s to study differentiation11. HL-60 undergoes cell cycle arrest and either myeloid or monocytic differentiation following stimulation; ATRA induces G1/G0-arrest and myeloid differentiation in HL-60 cells, while 1,25-dihydroxy vitamin D3 (D3) induces arrest and monocytic differentiation. Commitment to cell cycle arrest and differentiation requires approximately 48 hr of treatment, during which HL-60 cells undergo two division cycles.

Sustained mitogen-activated protein kinase (MAPK) pathway activation is a defining feature of ATRA-induced HL-60 differentiation. ATRA drives sustained activation of the Raf/MEK/ERK pathway, leading to arrest and differentiation12. Betraying a feedback loop, MEK inhibition results in the loss of ERK as well as Raf phosphorylation and the failure to arrest and differentiate in response to ATRA13. Retinoic acid (and its metabolites) are ligands for the hormone activated nuclear transcription factors retinoic acid receptor (RAR) and retinoid X receptor (RXR)14. RAR/RXR activation is necessary for ATRA-induced Raf phosphorylation13 in concert with the formation of an ATRA-induced signalsome complex at the membrane, which drives MAPK activation. While the makeup of the signalsome complex is not yet known, we do know that it is composed of Src family kinases Fgr and Lyn, PI3K, c-Cbl, Slp76, and KSR, plus transcription factors AhR and IRF115,16,17,18,19. Signalsome activity is driven by ATRA-induced expression of CD38 and putatively the heterotrimeric Gq protein-coupled receptor BLR120,21. BLR1 (also known as CXCR5), identified as an early ATRA (or D3)-inducible gene using differential display22, is necessary for MAPK activation and differentiation21, and drives signalsome activity. Studies of the BLR1 promoter identified a non-canonical RARE site consisting of a 17 bp GT box approximately 1 kb upstream of the transcriptional start that conferred ATRA responsiveness21. Members of the BLR1 transcriptional activator complex, e.g. NFATc3 and CREB, are phosphorylated by ERK, JNK or p38 MAPK family members suggesting positive feedback between the signalsome and MAPK activation23,24. BLR1 overexpression enhanced Raf phosphorylation and accelerated terminal differentiation, while Raf inhibition reduced BLR1 expression and ATRA-induced terminal differentiation25. In particular, Raf phosphorylation of the NFATc3 transcription factors at the BLR1 promoter enables transcriptional activation at the RARE by ATRA bound to RAR/RXR26. BLR1 knock-out cells failed to activate Raf or differentiate in the presence of ATRA25. Interestingly, both the knockdown or inhibition of Raf, also reduced BLR1 expression and functional differentiation25. Thus, the expression of signalsome components e.g., BLR1 was Raf dependent, while Raf activation depended upon the signalsome. A previous computational study of ATRA-induced differentiation of HL-60 cells suggested that the BLR1-MAPK positive feedback circuit was sufficient to explain ATRA-induced sustained MAPK activation, and the expression of a limited number of functional differentiation markers27. Model analysis also suggested that Raf was the most distinct of the MAPK proteins. However, this previous study developed and analyzed a complex model, thus leaving open the critical question of what is the minimal positive feedback circuit required to drive ATRA-induced differentiation.

In this study, we explored this question using a minimal mathematical model of the key architectural feature of ATRA induced differentiation of HL-60 cells, namely positive feedback between an ATRA-inducible signalsome complex and MAPK activation. The ATRA responsive signalsome-MAPK circuit was then used to drive a downstream gene expression program which encoded for the expression of intermediate and functional differentiation markers. The effective model used a novel framework which integrated logical rules with kinetic modeling to describe gene expression and protein regulation, while largely relying upon biophysical parameters from the literature. This formulation significantly reduced the size and complexity of the model compared to the previous study of Tasseff et al., while increasing the breadth of the biology described27. The effective model, despite its simplicity, captured key features of ATRA induced differentiation of HL-60 cells. Model analysis predicted the bistability of MAPK activation as a function of ATRA exposure; conformational experiments supported ATRA-induced bistability. Model simulations were also consistent with measurements of the influence of MAPK inhibitors, and the failure of BLR1 knockout cells to differentiate when exposed to ATRA. In addition, the expression of intermediate and phenotypic differentiation markers was also captured following ATRA exposure. Lastly, we showed through immunoprecipitation and inhibitor studies, that the guanine nucleotide exchange factor Vav1 is potentially a new ATRA-inducible member of the signalsome complex functioning as a regulator that contributes to signal amplification in the signalsome. Taken together, these findings when combined with other literature evidence, suggested that reinforcing feedback was central to differentiation programs generally, and necessary for ATRA-induced differentiation. The model answers a biologically important question that is not easily experimentally attacked, namely given the complexity of the signaling machine and the pathways it embodies, is there a critical small suite of molecules that are the action elements seminal to eliciting ATRA-induced cell differentiation and G0 arrest.

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).

Figure 1
figure 1

Schematic of the effective ATRA differentiation circuit. Above a critical threshold, ATRA activates an upstream Trigger, which induces signalsome complex formation. Signalsome activates the mitogen-activated protein kinase (MAPK) cascade which in turn drives the differentiation program and signalsome formation. Both Trigger and activated cRaf-pS621 drive a phenotype gene expression program responsible for differentiation. Trigger activates the expression of a series of transcription factors which in combination with cRaf-pS621 result in phenotypic change.

Full size image
Table 1 Myelomonocytc transcription factor connectivity used in the signal integration and phenotype modules.
Full size table
Table 2 Characteristic model parameters estimated from literature.
Full size table

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.

Figure 2
figure 2

Model analysis for ATRA-induced HL-60 differentiation. (A) BLR1 mRNA versus time following exposure to 1 \(\mu \)M ATRA at t = 0 hr. (B) cRaf-pS621 versus time following exposure to 1 μM ATRA at t = 0 hr. Points denote experimental measurements, solid lines denote the mean model performance. Shaded regions denote the 99% confidence interval calculated over the parameter ensemble. (C) Signalsome and cRaf-pS621 nullclines for ATRA below the critical threshold. The model had two stable steady states and a single unstable state in this regime. (D) Signalsome and cRaf-pS621 nullclines for ATRA above the critical threshold. In this regime, the model had only a single stable steady state. (E) Morphology of HL-60 as a function of ATRA concentration (t = 72 hr). Experimental data in panels A and B were reproduced from Wang and Yen25, data in panel E is reported in this study.

Full size image

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.

Figure 3
figure 3

Model simulation following exposure to 1 \(\mu \) M ATRA. (A) BLR1 mRNA versus time with and without MAPK inhibitor. (B) cRaf-pS621 versus time following pulsed exposure to 1 \(\mu \)M ATRA with and without BLR1. Solid lines denote the mean model performance, while shaded regions denote the 99% confidence interval calculated over the parameter ensemble. (C) Western blot analysis of phosphorylated ERK1/2 in ATRA washout experiments. Experimental data in panels A and B were reproduced from Wang and Yen25, data in panel C is reported in this study. The image of the raw gel for panel C is given in the Supplemental Materials.

Full size image

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 ATRA28,29,30,31. However, the bulk of the model parameters were taken from literature32 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.

Figure 4
figure 4

Model simulation of the HL-60 gene expression program following exposure to 1 μM ATRA at t = 0 hr. (A) Scaled CD38 and CD11b expression versus time following ATRA exposure at time t = 0 hr. (B) Scaled 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.

Full size image
Figure 5
figure 5

Model simulation of HL-60 cell-cycle arrest following exposure to 1 \(\mu \)M ATRA at t = 0 hr. (A) Predicted p21 and E2F expression levels for the best parameter set following ATRA exposure at time t = 0 hr. (B) Estimated fraction of HL-60 cells in G0 arrest following ATRA exposure at time t = 0 hr. Solid lines (or bars) denote the mean model performance, while the shaded region (or error bars) denotes the 95% confidence estimate of the polynomial model. Experimental data in panel B was reproduced from Jensen et al.31.

Full size image

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.

Figure 6
figure 6

Robustness of the HL-60 differentiation program following exposure to 1 \(\mu \)M ATRA at t = 0 hr. (A) 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.

Full size image
Figure 7
figure 7

Robustness of the HL-60 differentiation program following exposure to 1 \(\mu \)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.

Full size image

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 cells25,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 differentiation34 and their role in HL-60 differentiation has been investigated elsewhere15. 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. 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 cells17. 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.

Figure 8
figure 8

Investigation of a panel of possible Raf interaction partners in the presence and absence of ATRA. (A) Species identified to precipitate out with Raf: first column shows Western blot analysis on total Raf immunoprecipitation with and without 24 hr ATRA treatment and the second on total lysate. (B) The expression of species considered that did not precipitate out with Raf at levels detectable by Western blot analysis on total lysate. (C) Effect of the Raf inhibitor GW5074 on Raf interactions as determined by Western blot analysis of total Raf immunoprecipitation. The Authors note the the signal associated with Src was weak. (D) Cell Cycle distribution as determined by flow cytometry indicated arrest induced by ATRA, which was increased by the addition of GW5074. (E) Expression of the cell surface marker CD11b as determined by flow cytometry indicated increased expression induced by ATRA, which was enhanced by the addition of GW5074. (F) Inducible reactive oxygen species (ROS) as determined by DCF flow cytometry. The functional differentiation response of ATRA treated cells was mitigated by GW5074. GAPDH was used as a loading control. The black lines frame groupings from independent gels and each image is typical of three repeats. Images of the raw gels for panels A, B C and F are given in the Supplemental Materials.

Full size image

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 \(\pm \) 1%) compared with ATRA alone (61 \(\pm \) 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 arrest12. 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 differentiation35. 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 neutrophils20. 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 \(\pm \) 3.0% ROS induction (Fig. 8F). Cells treated with ATRA had a significantly increased TPA response, 53 \(\pm \) 7% ROS induction 48 hr post-treatment. Treatment with both ATRA and GW5074 statistically significantly reduced ROS induction (22 \(\pm \) 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 alterations36. 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 bistability37,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 progesterone39. 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α expression38. 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 bistable40. Still more biochemical networks important to human health, for example the human coagulation or complement cascades, also feature strong positive feedback elements41. 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 neutrophils42. PU.1, a member of the ets transcription factor family, is a well known regulator of granulocyte and monocyte development43. The relative level of PU.1 and C/EBPα is thought to control macrophage versus neutrophil cell fate decisions in granulocytic macrophage progenitor cells44. 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 neutrophils45. 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 infection42. 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 literature46,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-6033. 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 Raf48. 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 line49. CK2 can also phosphorylate Raf, although the literature has traditionally focused on S338 and not S621 or S25950. 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 activity51. Additionally, the association of Raf with 14-3-3 not only stabilized S621 phosphorylation, but also reversed the S621 phosphorylation from inhibitory to activating52. 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 cells53. 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 neutrophils54,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 induction55. 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 activation17. 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 STAT156,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 stimulus60. 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.

Figure 9
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 model17,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 \(\alpha \). 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 \(\alpha \) 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.

Full size image

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=\mathrm{1,2,}\ldots ,{\mathscr{G}}\), we modeled both the mRNA (\({m}_{j}\)), protein (\({p}_{j}\)) and signaling species abundance:

$$\frac{d{m}_{j}}{dt}={r}_{T,j}-(\mu +{\theta }_{m,j}){m}_{j}+{\lambda }_{j}$$
(1)
$$\frac{d{p}_{j}}{dt}={r}_{X,j}-(\mu +{\theta }_{p,j}){p}_{j}$$
(2)
$${\bf{g}}({p}_{1},\ldots ,{p}_{{\mathscr{G}}},\kappa )={\bf{0}}$$
(3)

where signaling species abundance was governed by the non-linear algebraic equations \({\bf{g}}({p}_{1},\ldots ,{p}_{{\mathscr{G}}},\kappa )={\bf{0}}\). The model parameter vector is denoted by \(\kappa \). The terms \({r}_{T,j}\) and \({r}_{X,j}\) denote the specific rates of transcription, and translation while the terms \({\theta }_{m,j}\) and \({\theta }_{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 \({\bar{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 \(0\le {u}_{j}\le 1\) depended upon the combination of factors which influenced the expression of gene \(j\). If the expression of gene \(j\) was influenced by \(\mathrm{1,}\ldots ,m\) factors, we modeled this relationship as \({u}_{j}={{ {\mathcal I} }}_{j}({f}_{1j}(\cdot ),\ldots ,{f}_{mj}(\cdot ))\) where \(0\le {f}_{ij}(\cdot )\le 1\) denotes a regulatory transfer function quantifying the influence of factor \(i\) on the expression of gene \(j\), and \({ {\mathcal I} }_{j}(\cdot )\) 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 expression62:

$${u}_{j}=\frac{{W}_{{R}_{1},j}+\sum _{n}\,{W}_{nj}{f}_{nj}}{1+{W}_{{R}_{1},j}+\sum _{d}\,{W}_{dj}{f}_{dj}}$$
(4)

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}_{1},j}\) denotes the weight of constitutive expression for gene \(j\). Given the formulation of the control law, the \({\lambda }_{j}\) term (which denotes the constitutive rate of expression of gene \(j\)) was given by:

$${\lambda }_{j}={\bar{r}}_{T,j}(\frac{{W}_{{R}_{1},j}}{1+{W}_{{R}_{1},j}})$$
(5)

The kinetic transcription term \({\bar{r}}_{T,j}\) was modeled as:

$${\bar{r}}_{T,j}={V}_{T}^{max}(\frac{{L}_{T,o}}{{L}_{T,j}})(\frac{{G}_{j}}{{K}_{T}+{G}_{j}})$$
(6)

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}\)), \({V}_{T}^{max}={k}_{T}({R}_{1})\). The \(({L}_{T,o}/{L}_{T,j})\) 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 \(({L}_{T,o}/{L}_{T,j})\) is a gene specific correction to the characteristic transcription rate \({V}_{T}^{max}\). If a gene expression process had no modifying factors, \({u}_{j}=1\). Lastly, the specific translation rate was modeled as:

$${r}_{X,j}={V}_{X}^{max}(\frac{{L}_{X,o}}{{L}_{X,j}})(\frac{{m}_{j}}{{K}_{X}+{m}_{j}})$$
(7)

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}_{X}^{max}={k}_{X}({R}_{2})\). 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.

Table 3 Sequence lengths from NCBI RefSeq database were used in the signal integration and phenotype modules128. The RNA sequence length used represents the total distance of transcription, and assume to be equal to the gene length.
Full size table

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\mathrm{,1}}\)) was estimated directly from the RAR/RXR abundance (\({x}_{u\mathrm{,1}}\)):

$${x}_{a\mathrm{,1}}\sim {x}_{u\mathrm{,1}}(\frac{\alpha \cdot {\rm{ATRA}}}{1+\alpha \cdot {\rm{ATRA}}})$$
(8)

where \(\alpha \) denotes a gain parameter; \(\alpha =0.0\) if ATRA is less than a threshold, and \(\alpha =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:

$$\frac{d{x}_{i}}{dt}={r}_{+,i}({\bf{x}},{\bf{k}})-(\mu +{k}_{d,i}){x}_{i}\quad \quad i=\mathrm{1,}\ldots ,{ {\mathcal M} }$$
(9)

The quantity \({x}_{i}\) denotes concentration of signaling species \(i\), while \( {\mathcal R} \) and \( {\mathcal M} \) denote the number of signaling reactions and signaling species in the model, respectively. The term \({r}_{+,i}({\bf{x}},{\bf{k}})\) denotes the rate of generation of activated species \(i\), while \(\mu \) 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:

$${x}_{i}^{\ast }\simeq \frac{{r}_{+,i}({\bf{x}},{\bf{k}})}{(\mu +{k}_{d,i})}\quad \quad i=\mathrm{1,}\ldots ,{ {\mathcal M} }$$
(10)

The generation rate was written as the product of a kinetic term (\({\bar{r}}_{+,i}\)) and a control term (\({v}_{i}\)). The control terms \(0\le {v}_{j}\le 1\) depended upon the combination of factors which influenced rate process \(j\). If rate \(j\) was influenced by \(\mathrm{1,}\ldots ,m\) factors, we modeled this relationship as \({v}_{j}={ {\mathcal I} }_{j}({f}_{1j}(\cdot ),\ldots ,{f}_{mj}(\cdot ))\) where \(0\le {f}_{ij}(\cdot )\le 1\) denotes a regulatory transfer function quantifying the influence of factor \(i\) on rate \(j\). The function \({ {\mathcal I} }_{j}(\cdot )\) is an integration rule which maps the output of regulatory transfer functions into a control variable. In this study, we used \({ {\mathcal I} }_{j}\in \{min,max\}\) and hill transfer functions63. If a process had no modifying factors, \({v}_{j}=1\). The kinetic rate of cRaf-pS621 generation \({\bar{r}}_{+,cRaf}\) was modeled as:

$${\bar{r}}_{+,cRaf}={k}_{+,cRaf}{x}_{s}(\frac{{x}_{cRaf}}{{K}_{+,cRaf}+{x}_{cRaf}})$$
(11)

where \({x}_{s}\) denotes the signalsome abundance, \({k}_{+,cRaf}\) denotes a characteristic activation rate constant, and \({K}_{+,cRaf}\) denotes a saturation constant governing cRaf-pS621 formation. In thus study, signalsome abundance was approximated by the abundance of the BLR1 protein; BLR1 expression is directly related to Raf nuclear translocation which in turn is related to activated signalsome. Thus, BLR1 is an indirect measure of the signalsome. The formation of cRaf-pS621 was regulated by only a single factor, the abundance of MAPK inhibitor, thus \({v}_{+,cRaf}\) took the form:

$${v}_{+,cRaf}=(1-\frac{{\rm{I}}}{{K}_{D}+{\rm{I}}})$$
(12)

where \({\rm{I}}\) denotes the abundance of the MAPK inhibitor, and \({K}_{D}\) denotes the inhibitor affinity.

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}_{ij}\in [\mathrm{0,100}]\); all cooperativity coefficients \({\eta }_{ij}\) appearing in the binding transfer functions \({f}_{dj}\) were bounded between \({\eta }_{ij}\in [\mathrm{0,4}]\); and all disassociation constants \({K}_{ij}\) appearing in the binding transfer functions \({f}_{dj}\) were bounded between \({K}_{ij}\in [\mathrm{0,1000}]\) (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\):

$${E}_{j}({\bf{k}})={({{ {\mathcal M} }}_{j}({t}_{-})-{\hat{y}}_{j}({t}_{-},{\bf{k}}))}^{2}+\sum _{i\mathrm{=1}}^{{{ {\mathcal I} }}_{j}}{({\hat{{ {\mathcal M} }}}_{ij}-{\hat{y}}_{ij}({\bf{k}}))}^{2}+\sum _{i\mathrm{=1}}^{{{ {\mathcal I} }}_{j}}{({{ {\mathcal M} }}_{ij}^{^{\prime} }-{y}_{ij}^{\text{'}}({\bf{k}}))}^{2}$$
(13)

The first term in Eq. (13) quantified the initial scale error, directly before the addition of ATRA. In this case, \({ {\mathcal M} }_{j}({t}_{-})\) (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 \({\hat{ {\mathcal M} }}_{ij}\) and \({\hat{y}}_{ij}\) 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 \({{\mathscr{T}}}_{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\le { {\mathcal M} }_{ij}^{\text{'}}\le 1\) was given by:

$${\hat{ {\mathcal M} }}_{ij}=({ {\mathcal M} }_{ij}-\mathop{{\rm{m}}{\rm{i}}{\rm{n}}}\limits_{i}{ {\mathcal M} }_{ij})/(\mathop{{\rm{m}}{\rm{a}}{\rm{x}}}\limits_{i}{ {\mathcal M} }_{ij}-\mathop{{\rm{m}}{\rm{i}}{\rm{n}}}\limits_{i}{ {\mathcal M} }_{ij})$$
(14)

where \({ {\mathcal M} }_{ij}^{\text{'}}=0\) and \({ {\mathcal M} }_{ij}^{\text{'}}=1\) describe the lowest (highest) intensity bands. A similar scaling was used for the simulation output. We minimized the total model residual \({\sum }_{j}{E}_{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 previously64.

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, \(\hat{{\mathscr{A}}}(t)\), as:

$$\hat{{\mathscr{A}}}(t)\simeq {a}_{0}+\sum _{i\mathrm{=1}}^{N-1}{a}_{i}{{\varphi }}_{i}({\bf{p}}(t),t)$$
(15)

where \({a}_{i}\) were unknown parameters, and \({{\varphi }}_{i}({\bf{p}}(t),t)\) 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:

$${{\varphi }}_{i}({\bf{p}}(t),t)={(\frac{t}{T}+\frac{p21}{E2F}{|}_{t})}^{(i-\mathrm{1)}}$$
(16)

The parameters \({a}_{0},\ldots ,{a}_{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% CO2 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 growth65. Experimental cultures were initiated at \(0.1\times {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 cRaf inhibitor GW5074 from Sigma-Aldrich (St. Louis, MO) was dissolved in DMSO with a stock concentration of 10 mM, and used at a final concentration of 2 μM. HL-60 cells were treated with 2 μM GW5074 with or without ATRA (1 μM) at 0 hr. This GW5074 dosage had a negligible effect on the cell cycle distribution, compared to ATRA treatment alone.

Immunoprecipitation and western blotting

For immunoprecipitation experiments, cells were lysed as previously described. 300 μg protein (in 300 μL total volume) per sample was pre-cleared with Protein A/G beads. The beads were pelleted and supernatant was incubated with Raf antibody (3 μL/sample) and beads overnight. All incubations included protease and phosphatase inhibitors in M-PER used for lysis with constant rotation at 4 °C. Bead/antibody/protein slurries were then washed and subjected to standard SDS-PAGE analysis as previously described17. All antibodies were purchased from Cell Signaling (Boston, MA) with the exception of α-p621 Raf which was purchased from Biosource/Invitrogen (Carlsbad, CA), and α-CK2 from BD Biosciences (San Jose, CA).

Morphology assessment

Untreated and ATRA-treated HL-60 cells were collected after 72 hr and cytocentrifuged for 3 min at 700 rpm onto glass slides. Slides were air-dried and stained with Wright’s stain. Slide images were captured at 40X (Leica DM LB 100 T microscope, Leica Microsystems).

Cell cycle quantification

1 × 106 cells were collected by centrifugation at 700 rpm and resuspended in 200 μl of cold propidium iodide (PI) hypotonic staining solution containing 50 μg/ml propidium iodine, 1 μl/ml Triton X-100, and 1 mg/ml sodium citrate (all Sigma, St. Louis, MO). Cells were incubated at room temperature for 1 h and analyzed by flow cytometry with a BD LSRII flow cytometer (BD Biosciences, San Jose, CA) using 488-nm excitation and emission collected through 550 long-pass dichroic and a 575/26 band-pass filters. Doublets were identified by a PI signal width versus area plot and excluded from the analysis.

CD11b quantification

Expression of cell surface differentiation markers was quantified by flow cytometry. 1 × 106 cells were collected from cultures and centrifuged at 700 rpm for 5 min. Cell pellets were resuspended in 200 μl 37 °C PBS containing 2.5 μl of phycoerythrin (PE)-conjugated CD11b (both from BD Biosciences, San Jose, CA). Following 1 h incubation at 37 °C cell surface expression levels were analyzed by flow cytometry. PE was excited at 488 nm and emission was collected through 505 long-pass dichroic and 530/30 band-pass filters.

Respiratory burst quantification (ROS)

1 × 106 cells were collected and centrifuged at 700 rpm for 5 min. Pellets were suspended in 500 μl of PBS containing 5-(and-6)-chloromethyl-2′,7′-dichlorodihydro-fluorescein diacetate acetyl ester (DCF, Molecular Probes, Eugene, OR) with either DMSO carrier blank solution or 12-0-tetradecanoyl-phorbol-13-acetate (TPA, Sigma, St. Louis, MO) suspended in DMSO. Cells were incubated for 25 min at 37 °C and then analyzed by flow cytometry as previously described. Oxidized DCF was excited by a 488 nm laser and emission was collected through a 505 long pass dichroic mirror and a 530/30 nm band pass filter. The shift in fluorescence intensity in response to TPA was used to determine the percent cells with the capability to generate inducible oxidative metabolites. Gates to determine percent positive cells were set to exclude 95 % of control cells not stimulated with TPA.

Availability of the model code

The signaling and gene expression model equations, and the parameter estimation procedure, were implemented in the Julia programming language. The model equations were solved using the ODE23s routine of the ODE package66. The model code and parameter ensemble is available under an MIT software license and can be downloaded from http://www.varnerlab.org.