Cell‐specific responses to the cytokine TGFβ are determined by variability in protein levels

Figures

• 

  Figure 1. Dynamics and variability of SMAD2 signaling in single cells AFluorescent reporter system to measure SMAD signaling dynamics in individual cells. SMAD2 was fused to the yellow fluorescent protein mVenus (YFP) under the control of the human ubiquitin C promoter (UbCp) with the selection marker G418 (Geneticin). As a nuclear marker, histone 2B (H2B) was fused to the cyan fluorescent protein mCerulean (CFP) under the control of UbCp with the selection marker hygromycin.BWestern blot analysis of endogenous and YFP‐tagged SMAD2 in a stable clonal reporter cell line and the corresponding parental cell line. Cells were stimulated with 100 pM TGFβ1 and analyzed after 3 h. GAPDH was used as a loading control.CWestern blot analysis of SMAD2 activation in SMAD2‐YFP reporter and parental MCF10A cells. Cells were stimulated with 100 pM TGFβ1, and SMAD2 phosphorylation was analyzed at indicated time points. GAPDH was used as a loading control.D, ELive‐cell time‐lapse microscopy images of MCF10A cells expressing SMAD2‐YFP following treatment with 100 pM TGFβ1 (D). White circles indicate the segmented nucleus, and the estimated cytoplasmic area is represented by red annuli. The indicated cell was…

• 

  Figure EV1. Dynamics and variability of SMAD2 signaling in single cells AWestern blot analysis of SMAD2 activation in SMAD2‐YFP reporter and parental MCF10A cells. Cells were stimulated with 100 pM TGFβ1 and SMAD2 phosphorylation was analyzed at indicated time points. GAPDH was used as a loading control. Independent experiments were quantified and normalized to maximum values. Error bars indicate standard deviation of biological repeats (n = 4). Note that phosphorylated YFP‐SMAD2 was at background levels at 0.25 h presumably due to lower expression levels.BExpression of SMAD target genes in parental and SMAD2 reporter cell lines. Expression kinetics of the SMAD target genes SMAD7, SnoN, and PAI‐1 upon 100 pM TGFβ stimulation were measured by qPCR in the indicated cell lines. β‐Actin was used as an internal control. Error bars indicate standard deviation of technical triplicates.CLive‐cell time‐lapse microscopy images of H2B‐CFP expression in MCF10A cells following treatment with 100 pM TGFβ. The same detail as in Fig 1D is shown. White circle indicates the segmented nucleus, and the estimated cytoplasmic area is represented by red annuli.D–FThe cell indicated in Fig 1D was tracked over 24 h.…

• 

  Figure 2. SMAD dynamics decompose into distinct signaling classes Time‐resolved analysis of SMAD2 nuclear to cytoplasmic localization for varying stimulus levels. Nuc/cyt SMAD2 ratios for eight individual cells (thin lines) as well as the population median (thick line) are shown. See Appendix Table S1 for number of cells analyzed.Median nuc/cyt SMAD2 ratio of cells stimulated with varying concentrations of TGFβ1 over 24 h. See Appendix Table S1 for number of cells analyzed.Individual cells were clustered into six signaling classes according to their time‐resolved nuc/cyt SMAD2 ratio using dynamic time warping (DTW). Each line represents the median over all cells of the indicated cluster. Cells stimulated with varying TGFβ1 concentrations as indicated in (B) were included in the analysis.Distributions of signaling classes depending on TGFβ dose.Silhouette plots of cells sorted according to TGFβ concentration (upper panel) or signaling classes (lower panel). Plots provide a graphical representation of how well the nuc/cyt SMAD2 ratios of each cell correspond to trajectories of other cells in its own group according to the cDTW measure. Positive silhouette scores indicate that SMAD2 responses…

• 

  Figure EV2. Clustering heterogeneous SMAD translocation dynamics using dynamic time warping A, BDose‐dependency of SMAD target gene expression. Parental MCF10A cells were stimulated with varying concentrations of TGFβ and SMAD7 (A) and PAI‐1 (B) expression was measured by qPCR at indicated time points. Error bars indicate standard deviation of technical triplicates.CIllustrative comparison of Euclidian distance and dynamic time warping (DTW). In the left panel, three single‐cell trajectories (red, blue, and black) with similar Euclidian distances are shown (see table). DTW performs a non‐linear alignment in time (middle panel) that compensates the temporal shift of the peaks in the red and blue trajectory (cyan lines). This leads to a lower DTW distance score compared to the DTW distance between the red and black or blue and black trajectories, which remain almost unchanged (see table).DDissimilarity matrix calculated pair‐wise by cDTW of single‐cell trajectories treated with varying TGFβ doses. Strength of the TGFβ stimulation increases from top left to bottom right.EHeat map of single‐cell time courses sorted according to hierarchical clustering. The corresponding dendrogram is shown on…

• 

  Figure EV3. Determining causes of heterogeneity in SMAD signaling Time‐resolved analysis of SMAD2 nuclear/cytoplasmic localization before and after stimulation with varying concentrations of TGFβ. Cells were imaged for 24 h before and 12 h after stimulation. Median nuc/cyt SMAD2 ratios for indicated stimulation levels are shown. See Appendix Table S1 for number of cells analyzed.Heat map of SMAD2 translocation in individual cells imaged for 24 h before stimulation with 25 pM or 5 pM TGFβ (similar to Fig 3A). Each horizontal line represents a single cell. Cells were sorted either by the time of the last division before stimulation (left) or by the amplitude of their response (right). Obviously, the time of last division and signaling responses are not correlated. See Appendix Table S1 for number of cells analyzed.Time of last cell division before stimulus for each signaling class after stimulation with the indicated concentrations of TGFβ. Single‐cell time courses were mapped onto the previously observed signaling classes (Fig 2C; see Appendix II.H). Distributions are overlapping; no significant trend in cell division time is observable. White lines indicate median; boxes include data between…

• 

  Figure 3. Heterogeneity in SMAD dynamics determined by cellular state Heat map of SMAD2 translocation in individual cells over time. Cells were imaged for 24 h before stimulation with 100 pM TGFβ1. Each horizontal line represents a single cell, and the nuc/cyt ratio is shown as indicated in the legend. Time of cell division is indicated by white marks. Cells were sorted either by the time of the last division before stimulation (left) or by the amplitude of their response (right). Cell cycle and response are not correlated. See Appendix Table S1 for number of cells analyzed.Mapping of SMAD2 translocation dynamics in individual cells to previously identified signaling classes (compare Fig 2C). Cells were imaged for 24 h before stimulation with varying TGFβ1 concentrations (Fig EV3A). For each trajectory, the most similar signaling class was determined using Euclidian distance to the median dynamics of the previously defined clusters (Fig 2C) as a similarity measure. Median nuc/cyt SMAD2 ratios for resulting mapped subpopulations are shown. See Appendix Table S1 for number of cells analyzed.Time of last cell division before stimulus for each signaling class (defined in B). Distributions are…

• 

  Figure 4. Mathematical modeling of TGFβ signaling AOutline of a tiered approach to model heterogeneous signaling in single cells (see text for details).BTopology of TGFβ pathway model. The oval shapes represent free receptors (blue), ligand (yellow), and ligand–receptor complex (gray). Extension “‐e” signifies endosomal species. Rectangles represent SMAD2 (blue), SMAD4 (green), and generic feedback regulator (yellow). Extensions “p” indicate phosphorylated and “n” nuclear species. Production and degradation are shown by phi symbols. State transitions and intercompartmental shuttling are indicated with arrows, enzyme catalysis with circle headed bars, and feedback inhibition with blunt headed bars.CCalibration of population‐average model by fitting to median SMAD2 translocation dynamics of cells stimulated with different TGFβ concentrations. Experimental data points correspond to Fig 2B. Model fits to other datasets are shown in Fig EV4 (see also Appendix Table S4); parameter values are provided in Appendix Table S5 and Table EV1.DMedium TGFβ degradation over time. Blue line shows the ligand concentration after an initial stimulus with 25 pM TGFβ1 as predicted by the best‐fit mathematical…

• 

  Figure EV4. Calibration of a population‐average model AModel calibration with median SMAD4 translocation dynamics of cells stimulated with different TGFβ concentrations. Solid lines represent model fit, filled circles measured data. Experimental data points correspond to Appendix Fig S2J.B, CDynamics of total cellular TGFβ receptor protein levels upon stimulation with 100 pM TGFβ. Dynamics of TGFβR1 (B) and TGFβR2 (C) are shown as fold change relative to unstimulated cells. Solid red lines represent model fit and circles data measured by Western blot analysis. For each observation point, data spread from triplicates is shown. Gray area indicates measurement noise estimated by the fitting algorithm.DRestimulation with 2.5 pM TGFβ. Stimulation was performed once at the beginning of the experiment (left), at the beginning and at 3 h (middle), or at the beginning and at 6 h (right) to replenish the ligand pool to its initial concentration. Solid lines represent model fits, and circles measured data points. Gray area indicates measurement noise estimated by the fitting algorithm.E, FModel calibration with median SMAD2 and SMAD4 translocation dynamics upon receptor inhibition. Cells were…

• 

  Figure 5. Modeling heterogeneous signaling dynamics in single cells AThe model of TGFβ signaling was fitted to six signaling classes observed upon stimulation with 100 pM TGFβ1. Median nuc/cyt SMAD2 ratios (circles) and model fits (solid lines) are shown.BFeatures of SMAD2 translocation dynamics. We considered the amplitude (E) and timing (T) of the first peak of nuclear translocation as well as the amplitude at 300 min (L) as a measure for the signaling activity upon adaptation of the pathway.C, DModel performance at varying noise levels. Heterogeneous signaling in response to a 100 pM TGFβ1 stimulus was simulated by signaling class‐based modeling (C) or a direct ensemble modeling (D) (see main text). Noise in protein expression is modeled as a combination of correlated and uncorrelated noise (see Appendix IV.B). The differences among single‐cell signaling features between model and data are calculated as sum of squared errors and normalized to the maximal deviation observed (color bar). For each combination of correlated and uncorrelated noise, 10,000 cells were simulated.EMeasured and predicted distributions of signaling features for two TGFβ stimuli (2.5 and 100 pM). A population of…

• 

  Figure EV5. Modeling heterogeneous signaling dynamics predicts negative feedback as major determinant of signaling classes Unimodal distribution of protein concentrations in artificial cell populations. Distributions of basal concentrations for TGFβR1, SMAD2, and SMAD4 are shown for a population of artificial cells assembled according to the observed proportion of signaling classes at 100 pM TGFβ using calibrated noise levels.Silhouette plots of artificial cells sorted according to TGFβ concentration (left panel) or mapped to experimentally observed signaling classes (right panel). Plots provide a graphical representation of how well the nuc/cyt SMAD2 ratios of each simulated cell corresponds to trajectories of other simulated cells in its own group. Positive silhouette scores indicate that SMAD2 responses are more similar to the own group, while negative scores signify that the corresponding trajectory is closer to any of the other groups. In general, signaling classes provide better separation than sorting according to stimulus levels.Transition between signaling classes depending on feedback strength. The response of a reassembled population of artificial cells to 5 and 25 pM TGFβ was…

• 

  Figure 6. Negative feedback determines decomposition into signaling classes Predicted distributions of signaling classes depending on TGFβ dose. Simulations were performed as described for Fig 5E. The simulated time courses were mapped onto the original clusters dynamics (Fig 2C) as described in Appendix II.H.Transition between signaling classes depending on stimulus strength. Same data as in (A). Black lines and their thickness indicate the direction and extent of transitions between signaling classes. Filled circle size indicates the proportion of artificial cells in the corresponding signaling class.Transition between signaling classes depending on feedback strength. The response of a reassembled population of artificial cells to 100 pM TGFβ1 was simulated with reduced feedback expression as indicated (see Appendix IV.D) and mapped to previously observed signaling classes. Black lines and their thickness indicate the direction and extent of transitions between signaling classes. Transitions with a probability below 1% were excluded for better visualization.Variation of model parameters across signaling classes. For 30 independent model fits to the experimentally observed signaling…