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 manipulated at a granular level, advanced therapies could be developed for a spectrum of disorders. In this study, we used computational and experimental tools to unravel the response of HL-60 human myeloblastic leukemia cells to Retinoic Acid (RA). HL-60 is an archetype model for the study of proliferation and differentiation mechanisms. HL-60 G0-arrests and differentiates when exposed to RA. By studying the differentiation and proliferation architecture of HL-60, which is limited to two possible differentiation paths, we will gain insight into the molecular programs of more complex stem cells.

A mathematical model of RA induced HL-60 cell-cycle arrest and differentiation was formulated and compared with protein- and Transcription Factor (TF)-array measurements. Previous studies have established that BLR1, a G-protein coupled receptor expressed following RA exposure, was required for RA-induced cell-cycle arrest and leads to atypical slow and persistent MAPK signaling. A draft HL-60 model was constructed and tested against BLR1 wild-type (wt), knock-out and knock-in HL-60 cell-lines with and without RA. The initial HL-60 model described the dynamics of 729 proteins interconnected by 1356 interactions using coupled Ordinary Differential Equations (ODEs). A modular System of Systems (SoS) strategy was used to model the HL-60 architecture. Separately identified MAPK, G1-cell-cycle, Calcium/G-protein and transcription/translation modules were integrated into a master two-compartment (cytosol and nuclear) HL-60 model. The submodels were each identified in different cell-lines. Ensembles of submodel parameters were estimated by comparing submodel simulations with in-vivo and in-vitro data sets. Missing connectivity between submodels was estimated from a manual review of >300 publications and from the STRING and NETWORKIN databases. The approximately identified modules were then assembled into the composite HL-60 network. HL-60 model parameters and initial conditions were estimated by comparing simulations with 26 in-vivo HL-60 data sets. Parameters from the submodel library were used as an initial guess for parameter fitting. An optimal step-size evolutionary algorithm was used to minimize the residual between simulations and the 26 in-vivo experimental constraints. A unique HL-60 parameter set was not identified. Rather, an ensemble of possible parameter sets was generated. Currently, the HL-60 parameter ensemble consists of 876 possible parameter sets. Of the 1356 parameters, 80% are constrained with a Coefficient of Variation (CV) of less than 100%.

Initial HL-60 model simulations revealed the functional connectivity between MAPK and BLR1. As a proof of concept, simulations of MAPK activation following RA treatment in BLR1 wt, knock-out and knock-in HL-60 cells were compared with the study of Wang and Yen (J. Biol. Chem., 283:4375 - 4386, 2008). The working hypothesis of the initial simulation studies was that the correct model structure alone, in the absence of fine parameter tuning, was sufficient to capture qualitative trends. Using the ensemble of parameter sets, we estimated the variation of key proteins in different HL-60 backgrounds with and without RA. Consistent with the Wang and Yen data and previous studies, the initial HL-60 model was able to show up-regulation of BLR1 expression following RA exposure. The model also correctly captured the sustained activation of MAPK following the addition of RA. Lastly, the dynamics of arrest and differentiation markers were correctly captured in the simulations. The initial simulation studies led to several testable linkages between BLR1 expression and MAPK activation. Two possible paths mediated by PKCα and PLCγ are being explored using Immunoprecipitation (IP) and siRNA knockdown studies. When taken together, the initial HL-60 simulations demonstrated that our subnetwork library in combination with monte-carlo parameter sampling was sufficient to explore possible interactions between BLR1 and MAPK activation in RA-induced HL-60 differentiation. Our results suggested that qualitative insight into the properties of complex networks could be gleaned from simulation despite parameter uncertainty. Lastly, initial simulation studies suggested experimentally testable connections between BLR1 expression and MAPK activation, in particular, the connection between PKCα/PLCγ and MAPK activation.

Extended Abstract Status: Not Uploaded