Modeling and analysis of retinoic acid induced differentiation of uncommitted precursor cells

Abstract

Manipulation of differentiation programs has therapeutic potential in a spectrum of human cancers and neurodegenerative disorders. In this study, we integrated computational and experimental methods to unravel the response of a lineage uncommitted precursor cell-line, HL-60, to Retinoic Acid (RA). HL-60 is a human myeloblastic leukemia cell-line used extensively to study human differentiation programs. Initially, we focused on the role of the BLR1 receptor in RA-induced differentiation and G1/0-arrest in HL-60. BLR1, a putative G protein-coupled receptor expressed following RA exposure, is required for RA-induced cell-cycle arrest and differentiation and causes persistent MAPK signaling. A mathematical model of RA-induced cell-cycle arrest and differentiation was formulated and tested against BLR1 wild-type (wt) knock-out and knock-in HL-60 cell-lines with and without RA. The current model described the dynamics of 729 proteins and protein complexes interconnected by 1356 interactions. An ensemble strategy was used to compensate for uncertain model parameters. The ensemble of HL-60 models recapitulated the positive feedback between BLR1 and MAPK signaling. The ensemble of models also correctly predicted Rb and p47phox regulation and the correlation between p21-CDK4-cyclin D formation and G1/0-arrest following exposure to RA. Finally, we investigated the robustness of the HL-60 network architecture to structural perturbations and generated experimentally testable hypotheses for future study. Taken together, the model presented here was a first step toward a systematic framework for analysis of programmed differentiation. These studies also demonstrated that mechanistic network modeling can help prioritize experimental directions by generating falsifiable hypotheses despite uncertainty.

Fig. 1

Overview of BLR1-MAPK positive feedback loop driving RA induced HL-60 arrest and differentiation. RA signals are intercepted by a family of RAR/RXR nuclear receptors which in turn drive the expression of genes with RARE promoter elements. One key RA-regulated protein is BLR1. BLR1, a putative G protein-coupled transmembrane surface receptor, drives an atypical sustained MAPK signal which in turn activates the expression of genes required for the execution of the cell-cycle arrest and differentiation programs. MAPK also activates factors in the BLR1 transcriptional activator complex resulting in positive feedback.

Fig. 2

Parameter identification strategy. (A) Multiple Monte-Carlo trajectories were used to randomly explore parameter space. The simulation likelihood was used to generate a family of parameter sets used in the simulation study. We generated N = 2377 possible parameter sets and selected the 100 sets with the highest likelihood for inclusion in the ensemble. (B) Coefficient of Variation (CV; standard deviation of a parameter relative to its mean value) for the parameter ensemble used in this study. A small CV suggested a parameter was tightly constrained by the training data used for model identification. Black circles represent the CV values for the full N = 100 sets used in the ensemble. The gray circles indicate the CV values for a sub-ensemble (N = 47) selected from the main ensemble and used in the robustness analysis study. CV values were sorted from lowest to highest relative to the full ensemble.

Fig. 3

Model simulations over the parameter ensemble captured the sustained activation of MAPK following RA exposure (1 μM) at time=1 h. Dashed lines denote the simulation mean. Shaded regions denote one ensemble standard deviation. (A) Experimental and simulated levels of BLR1 mRNA following RA exposure. (B) Time profile of phosphorylated RAF1 activation following RA exposure. (C) Simulated versus measured phosphorylated MEK activation following RA exposure. (D) Simulated versus measured phosphorylated ERK following RA exposure. Data was adapted from Wang and Yen.

Fig. 4

The model recapitulated RA-induced feedback between BLR1 expression and MAPK activation. (A) Simulated BLR1 expression normalized to wildtype (WT) with Raf inhibition (−, 50% decrease in Raf initial condition) and overexpression (+, 50% increase in Raf initial condition) 48 h after the addition of RA. (B) Simulated phosphorylated Raf levels normalized to wildtype (WT) with BLR1 knockout (KO, BLR1 gene initial condition set to zero) and over-expression (+, 50% increase in BLR1 gene initial condition) 12 h after the addition of RA. (bottom panel) Corresponding model training data adapted from Wang and Yen. First row: effect of Raf siRNA (left) and overexpression (right) on the expression of BLR1 (Northern). Second row: effect of BLR1 knockout and overexpression of the level of phosphorylated Raf (S621).

Fig. 5

Computationally predicted markers of RA-induced phenotypic shift. (A) Predicted p21-CDK4-cyclinD complex formation was consistent with the percentage of G1/0-arrested cells (insert). (B) Effect of RA on Rb expression. Rb transcript (top) remains constant while Rb protein (bottom) decreases. Rb transcript consistent with Northern analysis (top insert) while Rb protein levels were consistent with Western analysis (bottom insert). G1/0-arrest data was reproduced from Yen et al., Exp. Cell Res., 165: 193–151 1986. Rb Northern data was reproduced from Yen et al., Eur. J. Cell Biol., 65: 103–113 1994. Rb Western analysis was performed in our lab as described in the Experimental section. Dashed lines denote the simulation mean. Shaded regions denote one ensemble standard deviation.

Fig. 6

Robustness analysis. Each non-zero initial condition (conserved species) was removed, the model was run to approximate steady state and RA was added at time = 1 h. The area under the curve was calculated for each model species. (A) Qualitative coupling results. Removed species are along the x-axis from lowest to largest impact and observed model species are along the y-axis from least to most effected. Blue or red markers depict a statistical decrease or increase, respectively, in the area under the curve within a 90% confidence interval. (B) Coupling coefficients (area under the curve from the simulation with species removed over wild-type simulation) for three markers of differentiation: phosphorylated ERK, p47phox expression and p21-CDK4-cyclin D complex. Red circles indicate knock-downs which demonstrated a statistical decrease in all three markers: (left to right) RAR, RXR, BLR1, NFATc3, RNAp, eIF4E, Oct1, CREB, 40s and 60s ribosomes, met tRNA, EIF2, PABP, eIF4A, B, G, H, and eIF1, 1A, 3, 5, 5B. (C) Dendrogram of knockout species. The distance metric was the Euclidean norm and the linkage function was the inner square product (variance minimization algorithm). Each additional cluster is chosen to reduce the variance (y-axis). The color-threshold was chosen to be 200 which is 50% of the remaining variance after the initial division. General species and/or functions are indicated below each colored group. (D) Distinguishability is defined as the magnitude of the orthogonal components for all knockout species considered. Species are ordered from largest to smallest magnitudes. Red markers indicate species which are statistically significantly above 5. Specific species are identified as shown. Error bars show one standard deviation over the parameter ensemble.

Fig. 7

The HL-60 network architecture exhibits scale free properties.

Fig. 8

Correlation between parameter sets in the HL-60 ensemble. Regions of red indicate high correlation, while blue regions denote low correlation.