Understanding the molecular basis of stem cell differentiation programs is one of the grand unmet challenges facing modern cell-biology. If these programs could be rationally manipulated, advanced stem-cell therapies could be developed for the treatment of a spectrum of human diseases, including cancer. For example, differentiation induction cancer chemotherapy has been used successfully to treat a sub-type of leukemia, Acute Promyelocytic Leukemia (APL) which occurs primarily in young adults. Differentiation induction chemotherapy uses Retinoic Acid (RA), a form of vitamin A, to cause growth arrest and differentiation of leukemic cells, thereby obviating the high toxicity typical of conventional chemotherapy. However, leukemias other than APL as well as other cancers are not-responsive to RA treatment. In this study, we integrated molecular network modeling and analysis techniques with molecular cell-biology to determine the critical regulatory molecules that interfere with differentiation induction therapies for cancers other than APL. The experimental system being analyzed was the HL-60 human myeloid leukemia cell-line. HL-60, an archetype experimental model for the study of differentiation programs since the late 1970s, is RA sensitive however it has a genetic background similar to RA-insensitive cells. Thus, HL-60 offers an opportunity to understand how the RA-response program functions for cancers other than APL.

A mathematical model of the molecular programs governing RA induced cell-cycle arrest and differentiation was formulated from literature and compared with protein- and Transcription Factor (TF)-array measurements taken from HL-60 cells exposed to RA. RA exposure induced the expression of two receptors, CD38 and BLR1, also known as CXCR5, early in the terminal differentiation process. Both CD38 and BLR1 contribute to MAPK signaling and propel RA-induced myeloid differentiation. A draft HL-60 mathematical model was constructed and tested against BLR1 wild-type (wt), knock-out and knock-in HL-60 cell-lines in the presence and absence of Retinoic Acid (RA). The initial HL-60 model architecture described the dynamics of 88 proteins which translated, when reversible interactions were considered, into 563 proteins and protein complexes interconnected by 951 interactions. A modular strategy was employed for the construction of the model. Separately identified MAPK, G1-cell-cycle, Calcium/Gq-protein and translation initiation modules were integrated into a master two-compartment (cytosol and nuclear) HL-60 network model. The master network described the regulated and constitutive expression of 17 proteins (Oct1, CREB2, ETS, BRN, IRF,SRPK, pRb, RhoGDI, P47phox, CD45, EIF2AK, SIIIp15, AAF, Stat5b, XRE and SAP) following a two-operator site promoter template. Activated transcription factors were assumed to reversibly bind operator site 1 (O1) which then allowed the recruitment of RNAp to operator site 2 (O2). Transcription of BLR1 was modeling by explicitly describing the reversible assembly of the CREB2:Oct1:NFATc3 transcription complex, the dissociation of CREB2 to activate the complex at O1 and then recruitment of RNAp at O2. The total concentration of the remaining 61 proteins was fixed by the initial condition of the simulation. Missing connectivity between pathway modules was estimated from a manual review of >100 publications and from the STRING and NETWORKIN databases. Also, when not experimentally derived, NETWORKIN was used to estimate kinase-target relationships and likely phosphorylation sites. The 563 Ordinary Differential Equations (ODEs) of the master model were generated using the UNIVERSAL code generation tool along with code to perform monte-carlo sensitivity analysis on the model parameters. Mass action kinetics were used to describe all protein-protein and protein-DNA interactions. The current HL-60 model has 951 parameters distributed among three types, association (kon), dissociation (koff) and catalytic (kcat). UNIVERSAL assigned heuristic values to unknown parameters where the mean binding affinity was micromolar and the average kcat was uniformly distributed around 500 s^-1.

In the absence of parameter tuning, simulations of MAPK activation following RA treatment in BLR1 wt, knock-out and knock-in HL-60 cells were compared with the recent study of Wang and Yen (J. Biol. Chem., 283:4375 - 4386, 2008). The objective of these initial simulation studies was to test whether or not the proper model structure (biological connectivity) was in place to capture qualitative trends in the Wang and Yen data. Thus, the heuristic model parameters generated by UNIVERSAL were used in these initial studies (however the initial conditions were adjusted). A single parameter set was used in all simulations. Consistent with the Wang and Yen data, the initial HL-60 network model was able to show up-regulation of BLR1 following RA exposure along with generation of ERKp. The also model correctly captured the direction of shifts in expression of BLR1 following RA-exposure in Raf+/Raf- HL-60 lines and the influence of a BLR1-/BLR1+ background upon MAPK activation (ERKp and MEKp) following RA-exposure. However, while the qualitative trends were consistent with data on a relative time-scale, the absolute time-scales were not in agreement. Moreover, some trends were not captured, e.g., the decrease of BLR1 mRNA at long exposure times. Interestingly, the assembly of the modular model components led to several testable linkages between BLR1 expression and MAPK. The route that was active in the initial simulations was Gq-protein activation of PKC via PLCg and DAG. BLR1 expression was predicted to drive IP3 and DAG formation through coupling with Gq-protein and PLCg activation. In the model DAG was assumed to be linked with the activation of Protein Kinase C (PKC) by helping recruit it to the plasma membrane. PKC was then assumed, consistent with literature, to be able to phosphorylate both Ras and Raf leading to MAPK activation. In the absence of parameter tuning, the BLR1-PKC connection could explain the activation of ERK following RA exposure but not the dependence of BLR1 expression upon MAPK activity. A second connection between MAPK and BLR1 was discovered in the initial analysis of the model. ERKp was found to activate the NFATc3 protein, a component in the transcription factor complex driving BLR1 expression. This connection could explain the dependence of BLR1 expression on Raf. Initial sensitivity analysis results confirm the critical role of Raf phosphorylation.

When taken together, the HL-60 modeling and sensitivity analysis results suggested experimentally testable hypothesis on the interconnections between BLR1 and MAPK activation. Understanding these relationships could be a critical first step to deciphering why RA-treatments fail in cancers other than APL. More generally, these studies demonstrate that dynamic molecular models of protein-protein and protein-DNA networks can be a useful tool to study complex programs in cell-biology.