Deciphering the history of ERK activity from fixed-cell immunofluorescence measurements

Abstract

The RAS/ERK pathway plays a central role in diagnosis and therapy for many cancers. ERK activity is highly dynamic within individual cells and drives cell proliferation, metabolism, and other processes through effector proteins including c-Myc, c-Fos, Fra-1, and Egr-1. These proteins are sensitive to the dynamics of ERK activity, but it is not clear to what extent the pattern of ERK activity in an individual cell determines effector protein expression, or how much information about ERK dynamics is embedded in the pattern of effector expression. Here, we evaluate these relationships using live-cell biosensor measurements of ERK activity, multiplexed with immunofluorescence staining for downstream target proteins of the pathway. Combining these datasets with linear regression, machine learning, and differential equation models, we develop an interpretive framework for immunofluorescence data, wherein Fra-1 and pRb levels imply long-term activation of ERK signaling, while Egr-1 and c-Myc indicate more recent activation. Analysis of multiple cancer cell lines reveals a distorted relationship between ERK activity and cell state in malignant cells. We show that this framework can infer various classes of ERK dynamics from effector protein stains within a heterogeneous population, providing a basis for annotating ERK dynamics within fixed cells.

Fig. 1. ERK activity and target genes are dose-responsive to Epidermal Growth Factors.

a Schematic of the experimental method. Live cells were imaged in 96-well plates for 19 hours and immediately fixed following the end of the time-lapse measurements. Plates were subsequently stained for antibody-based measurements. b Condition average response measurements for live-cell ERK biosensor (EKAR signal, shown in arbitrary units) with increasing concentrations of EGF. Data are presented as the mean of each condition (n_replicates = 3 for all conditions). c Condition average response measurements depicted as a heatmap. Each row represents the average EKAR FRET measurement for a condition, indicated by the color scale (yellow, high ERK; blue, low ERK). EGF concentration is indicated by colored triangles as in Fig. 1b. MEKi = MEK inhibitor PD0325901 (100 nM) (see Supplemental Table 1 for n_eplicates for each condition). d Plots of single-cell EKAR signals (in arbitrary units) for five representative cells in each indicated condition. Bold lines indicate the mean of all cells in one well of the condition. Red lines indicate the time points where treatments were added. e MCF10A cells immuno-stained with cyclic immunofluorescence. Each row depicts the same group of cells. Scale bar = 100 μm. f Nuclear quantification of cyclic immunofluorescence measurements from listed EGF condition. Plots indicate the median (bar), 25th/75th percentiles (box), and range of the data (whiskers). The dashed line indicates the median of vehicle control (Imaging Media). Variance-corrected t-tests (two-sided) were conducted by comparing each EGF-treated condition to control. n_replicates = 3. P values, relative to the 0 EGF condition, are shown on each distribution. Throughout the study, n_replicates is used to denote independent experimental replicates performed on different days.

Fig. 2. ERK target gene expression moderately correlates with features of ERK dynamics.

a Single-cell heatmap for EKAR FRET measurements and corresponding ETG intensity, each row represents one cell (n_cells = 97,960, n_replicates = 3). ETG expression colored by intensity of immunofluorescence (IF) measurements (log₁₀). b Schematic of ERK dynamic features used for analysis. Frequency was calculated by estimating the mean normalized frequency of the power spectrum of the EKAR FRET measurement time series for each cell. c Pearson correlation (r) between each ERK feature and each cyclic IF measurement, where single- cell values were used. d Pearson correlation (r) between single-cell ETG measurements and the EKAR FRET measurement at each timepoint from the live-cell experiment. e Spatial heatmap of EKAR (left) and ETG (right) measurements from a single well (control condition). Heatmap is organized by proximity of cells to each other so that neighboring cells in the well are plotted closer to each other in the heatmap. ETG colormap indicates the relative log intensity of data within each column; outliers in pERK column skew colormap towards red. (black = NA). Magenta box indicates cells pictured in f. White arrows indicate additional cells that recently activated ERK, which resulted in higher Egr-1 expression (right). f Images corresponding to the cells plotted within the magenta box in e, representing an example of an association that was observed consistently across all 3 experimental replicates. All panels shown are registered images of the same cells, with the scale bar indicating 50 μm shown in the Hoechst-stained image.

Fig. 3. Convolutional neural network identifies non-linear signal transmission.

a For each ETG, three types of prediction models were separately trained. Top: Simplified schematic of convolutional neural network architecture containing two convolutional layers and three fully connected layers. Bottom left: Multiple variable regression where ERK activity at each time point is considered as a predictor variable (TS linear). Bottom right: Multiple variable linear regression where nine features of ERK activity are considered as predictor variables (Featurized linear). b Top: Bar plot indicating R² for three models used to predict ETG levels. Bottom: Bar plot indicating mean square error for three models used to predict ETG levels. Error bars represent standard error calculated using values from each fold of the 5-fold cross-validation partitions. c Scatter plot of the predicted and observed values of the CNN trained on all 190 timepoints (19 hr). The data represent standardized (z-scored) values. d Feature attribution heatmap showing the importance of each timepoint in the CNN model trained on 150 timepoints (15 h). Colormap represents relative values within each row. c and d represent the validation set of the first 5-fold partition.

Fig. 4. ERK target gene expression predicts history of ERK activation.

a Single-cell regression showing the coefficient of determination (R2) of linear regression models that use ETGs to predict each ERK feature. 10-fold cross-validation was conducted to retrieve the best test-set model. This model was then fit on the full dataset. “All” indicates multiple regression models using all ETGs as predictors. b Scatter plots of single-cell regression models showing line of best fit. Color indicates relative density of the data. c Scatter plot showing each cell’s predicted (x-axis) vs true (y-axis) value in the multiple linear regression (MLR) models. d Results of adding predictors to MLR models. Color of each point indicates which predictor was added at each step. e Average values were calculated for all cells with the same condition. These values were then used to fit regression models that predict each ERK feature using ETGs. f Scatter plots showing line of best fit and confidence intervals for condition average regression models. Each dot indicates the average of one condition.

Fig. 5. Cancer cell types display deficiencies in processing ERK dynamics.

a Condition average responses for live-cell ERK biosensor (EKAR) under four conditions. Data are presented as the mean of each condition (n_{well replicates} = 3). b Single-cell response plots to indicated condition. The bold line indicates the average of all cells in one well of the condition. c Pearson correlation (r) between single-cell protein measurements and the EKAR FRET measurement at each timepoint from the live-cell experiment. d Single-cell regression showing the coefficient of determination (R²) of linear regression models that use protein levels to predict each ERK feature. 10-fold cross-validation was conducted to retrieve the best test-set model. This model was then fit on the full dataset. “All” indicates multiple regression models using all proteins as predictors. e Percentage of cells classified as pRb-positive or pRb-negative; intensity thresholds identified individually for each cell line. Independent t-tests (two-sided) were conducted by comparing each condition to the Imaging Medium control. n_replicates = 3. p-vals (left to right): 0.000527, 0.000002, 0.033124, 0.020208, 0.001571. Error bars: standard deviation; a.u., arbitrary units.

Fig. 6. Inferring spatiotemporal ERK patterns using classification models.

a Average ERK activity in each cluster identified by k-means clustering of EKAR time series data in single cells. b Box plot showing median, 25th/75th quartiles, and range of ETG intensity in each cluster. One-way ANOVA test was conducted to compare the means of each group to each other. n_replicates = 3. *pval < 0.05, compared to indicated group(s). p-values indicated in Sup. Table 4. c AdaBoostM2 algorithm was trained to predict the cluster ID of each cell using its ETG measurements as predictors. The model with the best test-set performance, using 10-fold cross-validation, is shown. d Predictor importance estimates of each ETG in the model shown in c. e Comparison of model performance as a function of region size. Red: Single-cell models (as shown in a–d) for each cell line. Blue: Models for hexagonal regions (radius = 50 μm) of cells, where clustering and AdaBoostM2 models were performed on average ERK signaling and ETG expression values for each region. Green: Models generated using the entire image (702 μm by 785 μm), using a similar training procedure, except a standard decision tree was used due to the low sample size (n_samples: MCF10A₂₄₂, HCC827₇₁, A549₇₂, MCF7₆₈). Error bars: standard dev. of test-set accuracy across 10-fold cross-validation. Center of error bar: mean. MCF10A models were trained with 8 predictors; cancer cell models were trained with 19 predictors. f MCF10A images (Hoechst) overlayed with inferred signaling histories using single-cell (top) or regional (bottom) models. Left: Cells treated with EGF for 16 hours. Right: Cells treated with EGF for 12 hours, then treated with 100 nM MEKi for 4 hours. Some cells remain unlabeled due to incomplete predictor data. Legend: Bold lines indicate mean ERK activity in observed clusters (as in Fig. 6a); shaded regions indicate 25th and 75th percentiles.

Fig. 7. Mathematical model identifies limits of ERK activity prediction method.

a Ordinary differential equation model representing ERK-dependent modification of a transcription factor (TF), expression of mRNA, and expression of a protein (P) product. Superscript P denotes phosphorylation of a molecule. Lowercase k’s indicate rate parameters, uppercase K indicates a dissociation constant for feedback effects. Clock icons indicate a time delay (τ). b 10,000 cells were randomly picked from our main dataset. We simulated the gene response of 1000 genes for each cell using our experimentally collected EKAR measurements. c 1000 simulated genes (sim-ETGs) were randomly assigned the listed rate parameters while other rate parameters in the model remained constant. d Coefficient of determination (R²) of single variable linear regression models using each sim-ETG to predict the average ERK activity in each cell. e Pearson correlation (r) between each sim-ETG (rows) measurement and the EKAR values at each timepoint from the live-cell experiment. f Top: one ERK activity trace from one cell in the dataset. Rest: gene expression response of the top 5 predictors of the mean ERK activity. g Multiple regression models fit to predict each ERK feature. For each ERK feature prediction model, a single sim-ETG was added as a predictor at each step. To determine the order of sim-ETGs to add, we performed single regression and ranked sim-ETG by their ability to individually predict each feature.