Pathway crosstalk enables cells to interpret TGF-β duration

The detection and transmission of the temporal quality of intracellular and extracellular signals is an essential cellular mechanism. It remains largely unexplored how cells interpret the duration information of a stimulus. In this paper, we performed an integrated quantitative and computational analysis on TGF-β induced activation of SNAIL1, a key transcription factor that regulates several subsequent cell fate decisions such as apoptosis and epithelial-to-mesenchymal transition. We demonstrate that crosstalk among multiple TGF-β activated pathways forms a relay from SMAD to GLI1 that initializes and maintains SNAILl expression, respectively. SNAIL1 functions as a key integrator of information from TGF-β signaling distributed through upstream divergent pathways. The intertwined network serves as a temporal checkpoint, so that cells can generate a transient or sustained expression of SNAIL1 depending on TGF-β duration. Furthermore, we observed that TGF-β treatment leads to an unexpected accumulation of GSK3 molecules in an enzymatically active tyrosine phosphorylation form in Golgi apparatus and ER, followed by accumulation of GSK3 molecules in an enzymatically inhibitive serine phosphorylation in the nucleus. Subsequent model analysis and inhibition experiments revealed that the initial localized increase of GSK3 enzymatic activity couples to the positive feedback loop of the substrate Gli1 to form a network motif with multi-objective functions. That is, the motif is robust against stochastic fluctuations, and has a narrow distribution of response time that is insensitive to initial conditions. Specifically for TGF-β signaling, the motif ensures a smooth relay from SMAD to GLI1 on regulating SNAIL1 expression.


INTRODUCTION
Cells live in a state of constant environmental flux and must reliably receive, decode, integrate and transmit information from extracellular signals such that response is appropriate. [1][2][3][4] Dysregulation of signal transduction networks leads to inappropriate transmission of signaling information, which may ultimately lead to pathologies such as cancer. Therefore, a central problem in systems biology has been to untangle how quantitative information of cellular signals is encoded and decoded. In general cells respond to one or more properties of a stimulus, such as its identity, strength, rate of change, duration and even its temporal profile. [5][6][7][8][9][10][11] There are extensive studies on the dose-response curves to reveal how cells respond differentially to a signal with different strength. In comparison, how cells respond to the temporal code of signals is less studies, and its physiological relevance receives much attention recently since most extracellular signals exist only transiently and cellular responses show dependence on signal duration. [12][13][14][15][16] Transforming growth factor-β (TGF-β) is a secreted protein that regulates both transient and persistent cellular processes such as proliferation, morphogenesis, homeostasis, differentiation, and the epithelial-to-mesenchymal transition (EMT). [17][18][19][20][21] Because it plays essential roles in developmental and normal physiological process, and its dysregulation is related to cancer, fibrosis, inflammation, Alzheimer's disease and many other diseases, the TGF-β signaling pathway has been probed extensively as a potential pharmaceutical target. 22,23 Several quantitative studies have expanded our knowledge on how the TGF-β-SMAD signaling pathway transmits the duration and strength information of the signal. [24][25][26][27][28] TGF-β can activate both SMAD-dependent and multiple SMADindependent pathways, which then converge onto some downstream signaling elements. It is unclear how cells transmit and integrate information of the TGF-β signaling distributed among these parallel pathways. Addressing this question requires studies beyond the TGF-β/SMAD axis as in earlier work, where quantifying SMAD proteins serves as the fundamental readout. [24][25][26] Here, we focused on expression of SNAIL1, which is such a downstream target and plays a key role in regulating a number of subsequent processes. Our results confirmed that crosstalk between the SMAD-dependent and independent pathways is key for cells to decode and transmit temporal and contextual information from TGF-β. We posit that the mechanism may be a central mechanistic signal transduction process as many signaling pathways share the network structure.

RESULTS
Network analysis reveals a highly connected TGF-β signaling network Through integrating the existing literature, we reconstructed an intertwined TGF-β-SNAIL1 network formed with SMAD-dependent and SMAD-independent pathways ( Supplementary Fig. S1). For further studies we then identified a coarse-grained network composed of a list of key molecular species (Fig. 1, and Supporting text for details). Along the canonical SMAD pathway, TGF-β leads to phosphorylation of SMAD2 and/or SMAD3 (pSMAD2/3), followed by nuclear entry after recruiting SMAD4 and forming the complex. The complex acts as a direct transcription factor for multiple downstream genes including SNAIL1 and I-SMAD. 24,29 I-SMAD functions as an inhibitor of pSMAD2/3, thus closes a negative feedback loop. TGF-β also activates GLI1, a key component of the Hedgehog pathway, both through transcriptional activation by pSMAD2/3, and through suppressing the enzymatic activity of glycogen synthase kinase 3 (GSK3). The latter is constitutively active on facilitating GLI1 and SNAIL1 protein degradation in untreated epithelial cells, 30,31 thus suppressing GSK3 is expected to lead to GLI1 and SNAIL1 protein accumulation. Other non-SMAD signaling pathways, such as MAPK, ERK, et al. may also impact on SNAIL1 expression but to a less extent. 29, 32 We represented them as "others" in the model without further explicit treatment within the period of TGF-β treatment studied here. Therefore, the network integrates multiple feedforward loops that converge at the regulation of SNAIL1 transcription. In the following sections, we will examine the functional roles of individual pathways in the network using several human cell lines.
The canonical TGF-β/SMAD pathway initializes a transient wave of SNAIL1 expression First we examined the TGF-β/SMAD/SNAIL1 pathway (Fig. 2a) by treating human MCF10A cells with recombinant human TGF-β1, and performing multicolor immunofluorescence (IF) using antibodies directed against pSMAD2/3, SNAIL1. As expected from the pSMAD/I-SMAD negative feedback loop, pSMAD2/3 proteins accumulated in the nucleus transiently, peaking at around 12 h after TGF-β1 treatment, followed by a decrease by 24 h (Fig. 2b, c). We confirmed the transient pSMAD2/3 dynamics by sampling 1100-2600 cells at each time point (Fig. 2c). The result is also consistent with reports in the literature. 24,26,33 Nuclear SNAIL1 concentration rose concurrently with pSMAD2/3 (Fig. 2b, c), then there was a transient dip at 24 h, followed by another increase then a persistant elevation for one week. 34 Next, we investigated the function of phosphorylated SMAD2/3 on promoting snail1 transcription during TGF-β treatment. In addition to adding TGF-β, we treated MCF10A cells with an inhibitor LY2109761, which prevents SMAD2/3 phosphorylation through inhibiting TGF-β receptor kinase activity (Fig. 2d). Without the inhibitor, the SNAIL1 mRNA showed the two-wave dynamics consistent to that of the protein. When the inhibitor was added concurrently with TGF-β treatment, the SNAIL1 mRNA level was reduced to~9% of that of the control experiment (without inhibitor) by day 3. This result is consistent with previous observation that directly blocking SMAD2/3 phosphorylation or pSMAD activation at the early stage of TGF-β treatment depletes snail1 expression significantly (1)(2)(3)(4), and indicates that indeed pSMAD2/3 are required for SNAIL1 initial activation. However, the SNAIL1 mRNA level remained~70% when the inhibitor was added 48 h after initiation of TGF-β treatment (when nuclear pSMAD2/3 concentration has dropped to a minimum). Furthermore we constructed a mathematical model that contains only the TGF-β/ SMAD/SNAIL1 pathway, and performed a thorough parameter space search using a multi-configuration Monte Carlo algorithm ( Supplementary Fig. S2). The search revealed regions of the parameter space that quantitatively reproduced the transient pSMAD2/3 dynamics, but not the two-wave dynamics of SNAIL1 expression (Fig. 2e). This computational result further confirmed that pSMAD2/3 is less essential for the second wave of SNAIL1.
Furthermore, this SNAIL1 dynamics is not cell type specific as equivalent two-wave dynamics were seen for SNAIL1 mRNA in MCF7 and A549 cells (Fig. 2f). Similar to that of MCF10A, it is more effective on inhibiting SNAIL1 mRNA by adding LY2109761 together with TGF-β than later (Fig. 2g). The impact of SMAD phosphorylation inhibitor on A549 is less than that on MCF10A or MCF7 at either early or late inhibition, which could be due to the higher level of EMT-related factors in A549. 35 In total, these results reveal that pSMAD2/3 is essential for the early phase of SNAIL1 activation, but is less important for the secondary phase elevation and persistence of SNAIL1 expression/localization. GLI1 contributes to activating the second wave of SNAIL1 The regulatory network suggests that GLI1 may be responsible for the second wave of SNAIL1 (Fig. 3a). To test this hypothesis, we performed microscopy studies of SNAIL1-GLI1 using MCF10A cells. The distribution of SNAIL1 found in this study (Supplementary Fig. S3a) was consistent with those from the pSMAD2/3-SNAIL1 studies. Elevated and sustained expression of GLI1 under TGF-β treatment (Fig. 3b, c) was clearly evident. More interestingly GLI1 also showed an unexpected multi-phasic dynamic. Around 8 h after TGF-β treatment, cytosolic GLI1 concentration started to increase. At 12 h when SMAD activities decreased toward basal levels there was a clear accumulation of GLI1 in the nucleus, which continued to increase through day 2. Notably, at this time point cells expressing a high level of nuclear SNAIL1 consistently showed high nuclear GLI1 concentrations ( Supplementary Fig.  S3a). Expanding the mathematical model of the network to Fig. 2a also reproduced the temporal dynamics of pSAMD2/3 and SNAIL1 ( Supplementary Fig. S3b), supporting the role of GLI1 as the activator of the second wave of SNAIL1.
If GLI1 is mainly involved only in the later maintenance of SNAIL1 expression, it is reasonable to predict that inhibiting GLI1 activity, either at the onset of or at some subsequent time after TGF-β treatment, would have minimal effect on the pSMAD2/3 induced initial wave of SNAIL1 expression. However, GLI1 inhibition would severely reduce the second wave of SNAIL1 expression. Indeed this was what observed experimentally. When GLI1 inhibitor GANT61 was added together with TGF-β at the beginning of the experiment, the SNAIL1 mRNA level was reduced to be 55% (at 12 and 24 h), 12% (at 48 h) and 7% (at 72 h) compared to that without inhibition at the corresponding time points (Fig. 3d). In another experiment adding the inhibitor 48 h after TGF-β treatment also reduced the mRNA level measured at 72 h to be 25% (Fig. 3e). These results are qualitatively different from those with the SMAD inhibitor (Fig. 2d).
To confirm that GLI1 activation is not restricted to the MCF10A cell line, we also examined MCF7 and A549 cells with TGF-β treatment. We observed similar increased and sustained GLI1 expression, albeit with initial slight downregulation before 12 h, possibly due to cell line specific activation of some GLI1 inhibition pathways (Fig. 3f). Furthermore, early and late GLI1 inhibition lead to a reduction of the SNAIL1 mRNA level to be 13 and 22% for MCF7 cells, and to a less extent of 57 and 66% for A549 cells, respectively (Fig. 3g). Additionally, increased GLI1 expression after TGF-β treatment has been found for multiple liver cancer cell lines. 36 In toto these results support the role of GLI1 as a signaling relay from pSMAD2/3 to SNAIL1. Reconstructed literature-based pathway crosstalk for TGF-β induced SNAIL1 expression. The node "others" refer remaining SNAIL1 activation pathways that have minor contributions to the time window under study and thus are not explicitly treated GSK3 in a phosphorylation form with augmented enzymatic activity accumulates at endoplasmic reticulum and Golgi apparatus Next, we hypothesized that GSK3 is fundamental to the observed multi-phasic GLI1 dynamic (Fig. 3b). Most published studies suggest that GSK3 is constitutively active in untreated cells, facilitating degradation of SNAIL1 and GLI1; TGF-β treatment leads to GSK3 phosphorylation and inactivation, which leads to an accumulation of SNAIL1 and GLI1. 37,38 Initially we tested whether the above mechanism is sufficient to explain the multi-phasic GLI1 dynamics. We treated MCF10A cells in the absence of TGF-β with a GSK3 activity inhibitor. Given the Fig. 2 The SMAD proteins induce the first wave of SNAIL1. a Canonical SMAD-dependent pathway for TGF-β activation of SNAIL1 highlighted from the network in Fig. 1. b Two-color immunofluorescence (IF) images of pSMAD2/3 and SNAIL1 of MCF10A cells induced by 4 ng/ml TGF-β1 at various time points. The scale bar is 10 μm and is the same for other IF images in this paper. c Distributions of nuclear pSMAD2/3 and SNAIL1 concentrations quantified from the IF images. Red vertical lines indicate the mean value of the distributions at time 0, and blue vertical lines represent that at 12 h (for pSMAD2/3) or at 48 h (for SNAIL1), respectively. The number marked in each figure panel is the number of randomly selected cells used for the analysis. Throughout the paper we report fold changes of concentration and amount relative to the mean basal value of the corresponding quantity. d Effects of early (added together with TGF-β) and late (48 h after adding TGF-β) pSMAD inhibition on the SNAIL1 mRNA level in MCF10A cells. e Thorough parameter space search confirmed that with the model in panel a one can fit the pSMAD2/3 dynamics, but not the two-wave SNAIL1 dynamics. The experimental data are shown as violin plots with the medians given by black bars. Solid curves are computational results with parameter sets sampled from the Monte Carlo search, and the red curves are the bestfit results. f Fold change of SNAIL1 mRNA levels in MCF7 and A549 cells measured with quantitative RT-PCR after TGF-β1 treatment. g Fold change of SNAIL1 mRNA levels measured with quantitative RT-PCR at 72 h after TGF-β1 (T) treatment. For early inhibition (T + I) the inhibitor was added at the time of starting TGF-β1 treatment. For late inhibition (T−/+I) the inhibitor was added 48 h (for MCF7) and 24 h (for A549) after starting TGF-β1 treatment, respectively. The inhibition results were compared to the TGF-β treatment (T) result at the same time point above mechanism, one should expect the GSK3 inhibitor to promote both GLI1 and SNAIL1. In our experiment, SNAIL1 did increase in the nucleus and even more in the cytoplasm due to inhibition of GSK3-dependent SNAIL1 degradation, but there was no noticeable change in GLI1 expression in either nucleus or cytoplasm (Fig. 4a), suggesting additional signaling mechanisms may be involved.
Besides the inhibitory serine phosphorylation (S21 in GSK-3α and S9 in GSK-3β), previous studies showed that tyrosine (Y279 in GSK-3α and Y216 in GSK-3β) phosphorylation leads to augmented enzymatic activity of GSK3. 39 As a convenience when discussing the three forms of GSK3, we refer the enzymatically active unphosphorylated form and the more active tyrosine phosphorylated form as "GSK3 A " and "GSK3 AA ", respectively, and the inactive serine phosphorylation form as "GSK3 D ". Also we reserve "GSK3" for the total GSK3. As expected, microscopy studies showed an increased concentration of GSK3 D peaking around 12 h after TGFβ treatment (Fig. 4b, Supplementary Fig. S4a). Large cell-to-cell variations in the concentration of GSK3 D were observed, however, the abundance of cytosolic and nuclear GSK3 D were essentially equivalent (the expression ratio was close to one) for cells without TGF-β treatment ( Supplementary Fig. S4b). This observation corroborates earlier report that the serine phosphorylation does not affect GSK3 nuclear location. 40 TGF-β treatment led to transient deviation of this ratio from equivalence, reflecting additional active and dynamic regulation of GSK3 including covalent modification, location and protein stability. Specifically prior to inhibitory serine phosphorylation we observed transient GSK3 AA accumulation in the perinuclear region peaking at 8 h (Fig.  4b, Supplementary Fig. S4a). Close examination of higher magnification confocal images revealed that the GSK3 AA formed clusters in the endoplasmic reticulum (ER) and Golgi apparatus, but not associated with actin filaments (Fig. 4c, Supplementary Movie 1 & 2). Given that a function of active GSK3 is to modify target proteins post-translationally, our observation suggests an unreported role for GSK3 AA accumulating at the ER and Golgi apparatus as to modify newly synthesized proteins before their release to the cytosol. Specifically previous studies showed that in mammalian cells a scaffold protein SUFU binds to GLI to form an inhibitory complex; SUFU phosphorylation by GSK3β prevents the complex formation, and exposes the GLI1 nuclear localization sequence. 41 This mechanism explains the observed increase of free GLI1 in the cytosol followed by nuclear translocation (Fig. 3b). Since the two phosphorylation forms, GSK3 AA and GSK3 D , coexist within single cells at defined time points, we performed co-immunoprecipitation and found that the probability of having the two GSK3 phosphorylation forms in one molecule was rare ( Supplementary Fig. S4c).
Contrary to our observation that TGF-β regulates GSK3 AA dynamics, other studies posit that GSK3 AA is not regulated by external cues. 42 To resolve this paradox, we measured the relative amount of different GSK3 forms through silver staining (Fig. S4d). Among the three forms, the overall percentage of GSK3 D increased from a basal level of 37-65% at 12 h after TGF-β treatment. In contrast, only a small fraction of GSK3 molecules assumed the GSK3 AA form and its overall abundance was stable over time (from~10% basal level to~13% at 8 h then back tõ 10% at 12 h after TGF-β treatment). Essentially GSK3 AA did not change in abundance but did change in localizations (homing to the ER and Golgi apparatus) to form a high local concentration, which imbue an important role in TGF-β signaling.
A temporal and compartment switch from active to inhibitory GSK3 phosphorylation smoothens the SMAD-GLI1 relay and reduces cell-to-cell heterogeneity on GLI1 activation Based on the above results, we constructed an expanded network for TGF-β induced SNAIL1 expression (Fig. 5a), which integrates a role for GSK3 and its temporal change of enzymatic activities in the cytosol and nucleus (Supplementary Fig. S4e). The model reproduces the multiphasic dynamics of GLI1 as well as that of pSMAD2/3 and SNAIL1 ( Supplementary Fig. S5a).
To understand the function of the early nuclear accumulation of GLI1 induced by GSK3 AA , it is important to recognize that GLI1 has a positive-feedback loop, and this network motif (Fig. 5b, left panel without the part in green) has characteristic sigmoidal shaped temporal dynamics, with the substrate concentration increasing slowly at first then accelerating with time until it approaches saturation (Fig. 5b, right panel, red curve). The response time, t R , defined as the time taken to reach a target concentration value (X) R , is highly sensitive to initial substrate concentration (X) 0 : in fact a slight increase in the initial concentration, Δ(X), can significantly shorten the response time (Fig. 5b, right panel, blue curve). In contrast, one can accelerate the response time with an expanded network (Fig. 5b including the green part) that the signal triggers a fast conversion of the substrate from a preformed inhibitory form (X I ) to active form X, effectively a boost of (X) 0 to (X) 0 + Δ(X). For a fixed Δ(X) a greater acceleration is seen in cells with lower initial concentrations (Fig. 5b, right panel, inlet figure). Consequently despite variations of their initial concentration (X) 0 , most cells within a population can reach (X) R by a targeted time point t T in a series of temporally regulated events such as cell differentiation and immune response. Indeed, many examples of this modified feedback loop motif exist. Figure S5b gives some examples involving members of intrinsically disordered proteins and inhibitors of DNA binding proteins, β-catenin and the STING motif for immune responses. In the present scenario the accelerated GLI1 dynamic ensures sufficient accumulation of GLI1 before nuclear pSMAD2/3 level decreases, essentially analogous to a relay race when the first runner can only release the baton after the second runner has grabbed it. Later when the GLI1 and SNAIL1 concentrations start to increase, the GSK3 A →GSK3 D conversion became necessary to reduce the rates of their degradation catalyzed by active GSK3. Interestingly, this conversion takes place concurrently with maximal concentration of nuclear pSMAD2/3, which activates GLI1 and SNAIL1 transcription. Furthermore, the small initial concentration boost does not affect another major function of the positive feedback loop, which is to robustly buffer temporal and strength fluctuations of signals ( Supplementary Fig. S5c). 43 To test the functional roles of GSK3 suggested above, we performed a series of GSK3 activity inhibition experiments. First, we pretreated MCF10A cells with GSK3 inhibitor SB216763, washed out the inhibitor then added TGF-β1 ( Supplementary  Fig. S5d). We predicted that the treatment would slow down GLI1 nuclear accumulation, and at later times decrease the overall increase of GLI1 and SNAIL1 compared to cells without GSK3 inhibitor. Indeed this was observed (Fig. 5c, TGF-β+/− GSK3_I). More interestingly, the scatter plots (Fig. 5d) show the distributions with and without the inhibitor are similar in cells with high GLI1, but in the presence of the inhibitor there is a population of non-responsive cells with low GLI1 and SNAIL1. This observation is consistent with model predictions that the GSK3-induced boost of initial GLI1 concentration leads to acceleration in the GLI1 and SNAIL1 dynamics, and this boost is more evident for cells with lower level of initial nuclear GLI1 (Fig. 5e). In a separate experiment ( Supplementary Fig. S5e), we did not wash out GSK3 inhibitor while adding TGF-β. In this case the inhibitor had opposite effects on GLI1 and SNAIL1 protein concentrations: it slowed down the initial release and translocation of GLI1 needed to accelerate the GLI1 accumulation, but also decreased GLI1 and SNAIL1 degradation that becomes pre-eminent when the proteins were present at high levels. Compared to the samples grown in the absence of the GSK3 inhibitor, we also observed slower and more scattered GLI1 nuclear accumulation and SNAIL1 increase on day 2, but by day 3 the overall levels of GLI1 and SNAIL1 were actually higher than the case without the inhibitor (Fig. 5d, TGF-β + GSK3_I).
The SMAD-GLI1 relay increases the network information capacity and leads to differential response to TGF-β duration Our results show that TGF-β1 signaling is effected through pSMAD2/3 directly with fast pulsed dynamics concurrently with a relay through GLI1 which has a much slower dynamics. The signaling ported by these two channels converges on SNAIL1 with a resultant two-wave expression pattern. To further dissect the potential functional interactions between these two pathways, we performed mathematical modeling and predicted that the two distinct dynamics allows cells to respond to TGF-β differentially depending on stimulus duration (Fig. 6a). Short pulses of TGF-β only activate pSMAD2/3 and the first wave of transient SNAIL1 expression. When the signal duration is longer than a defined threshold value, activation of GLI1 will lead to the observed second wave of SNAIL1 expression. We confirmed the predictions with MCF10A cells (Fig. 6b). Both TGF-β1 pulses with duration of 2 and 8 h activated pSMAD2/3 and the first wave of SNAIL1 expressions. However, only the 8-h but not the 2-h pulse activated sustained GLI1 and the second wave of SNAIL1 expression, similar to those with continuous TGF-β1 treatment.
Clearly cellular responses have different temporal profiles depending on the TGF-β duration, and one can use the information theory to quantify their information content. 9,10 In this study we utilized a more intuitive understanding of network function from an information encoding viewpoint. Consider the pSMAD complex, which has three coarse-grained states, high (H), medium (M), and low (L), and each of GLI1 and SNAIL1 has two states, H and L (Fig. 6c) Further modeling suggests that components in the network function cooperatively to encode the TGF-β information (Supplementary Fig. S6b). Increasing or decreasing the nuclear GSK3 enzymatic activity tunes the system to generate the second SNAIL1 wave with a higher or lower threshold of TGF-β duration, respectively, while changing the cytosol GSK3 enzymatic activity has the opposite effect. Upregulation of GLI1, or downregulation of I-SMAD, both of which have been observed in various cancer cells, also decrease the threshold for generating the second SNAIL1 wave. Therefore cells of different types can share the same network structure, but fine-tune their context-dependent responses by varying some dynamic parameters, and for a specific type of cells dysregulation of any of the signaling network components may lead to misinterpretation of the quantitative information of TGF-β signal.
We have shown that the SNAIL1 dynamics is TGF-β1 duration dependent. To further confirm that cells respond differentially to TGF-β1 with different duration, we measured the mRNA levels of another four genes, all of which respond to TGF-β1 (Supplementary Fig. S6c). 44,45 Gene FN1 codes for the cell motility related protein fibronectin. Its expression is activated even by the 2-h TGF-β1 pulse, and increases with longer TGF-β1. Gene CTGF, whose product is an extracellular matrix protein and related to cell motility, is activated at similar extent by both 2 and 8-h TGF-β1 pulses, and its expression level increases by additional 14-folds with continuous TGF-β1 treatment. Expressions of genes MMP2 and CLDN4, coding proteins related to mesenchymal extracellular hallmark and cell migration, increase only slightly (less than two folds) with either 2 or 8 h TGF-β1 pulse, compared to the more significant change under continuous TGF-β1 treatment. Therefore, Fig. 6 The TGF-β-SNAIL1 network permits detection of TGF-β duration and differential responses. a Model predictions that the network generates one or two waves of SNAIL1 depending on TGF-β duration. The red line overlaid on the heatmap is a sampling time of the shorttime TGF-β induction. The green line represents the long-time TGF-β treatment. b Single cell protein concentrations quantified from IF images of cells under pulsed and continuous TGF-β treatments. The solid lines divide the space into coarse-grained states with respect to the corresponding mean values without TGF-β treatments (=1). c Schematics of how cells encode information of TGF-β duration through a temporally ordered state space these downstream genes also show differential expression patterns depending on TGF-β1 duration, and cells activate different response programs correspondingly.

DISCUSSION
TGF-β is a multifunctional cytokine that can induce a plethora of different and mutually exclusive cellular responses. A significant open question is how cells interpret various features of the signal and make the cell fate decision. TGF-β can activate a number of pathways interconnected with multiple crosstalk points. Our studies reveal that this interconnection is essential such that components of the network can function coordinately and appropriately to interpret the temporal (time and duration) information from TGF-β.
pSMADs are major inducers for the first wave of SNAIL1 expression The two-wave dynamic of TGF-β-induced SNAIL1 expression has been observed in several cellular systems, 46,47 supporting the underlying relay mechanism discovered in this work. The first wave is fundamentally induced by pSMAD2/3, as evidenced from our SMAD inhibition experiments, and similarity between the dynamics of pSMAD2/3 and the first wave of SNAIL1. SNAIL1 may act as cofactor of pSMADs to induce other early response genes. 48 At later times the nuclear concentrations of pSMAD2/3 decrease though continue to contribute to SNAIL1 activation at a lower level.
GLI1 is a signaling hub for multiple pathways and temporal checkpoint for activating second-wave of sustained SNAIL1 expression GLI protein has been traditionally attributed to the canonical Hedgehog pathway. Here, we show that TGF-β induction of GLI1 relays the signal to induce SNAIL1. Consistent with the present study, Dennler et al. reported SMAD3-dependent induction of the GLI family by TGF-β both in multiple cultured cell lines, and in transgenic mice overexpressing TGF-β 30 Many other signals such as PGF, EGF can also activate the GLI family, and a GLI code has been proposed to integrate input from different pathways and lead to context-dependent differential responses. 49 Our results confirmed this role of GLI1 as an intermediate information integrator and transmitter, and suggest that TGF-β must act above a threshold value of duration to activate the second wave of SNAIL1. This temporal checkpoint prevents spurious SNAIL1 activation and subsequent major cellular fate changes.
GSK3 fine-tunes the threshold of the GLI1 checkpoint and synchronizes responses of a population of cells The functional switch from pSMAD2/3 to GLI1 relays information from TGF-β signaling beyond the initial induction of SNAIL1, and this relay is facilitated by a second relay from the active to the inactive phosphorylation form of GSK3 proteins. Active regulation of the abundance and nuclear location of GSK3 AA form has been observed in neurons. 50 In contrast to these earlier reports we observed an accumulation of GSK3 AA in the ER and Golgi apparatus. Mechanistically this may be caused by redistribution of cytosolic GSK3 AA , or a simple accumulation of de novo synthesized and phosphorylated GSK3 proteins. The overall consequence is an increase in local GSK3 enzymatic activity, which forms part of the GSK3 switch that smooths the pSMAD2/3-GLI1 transition and the duration threshold of TGF-β pulse that generates the second wave of SNAIL1.
This seemingly simple process, which accelerates the response time through transient and minor increases in the initial concentration of a molecular species subject to positive feedback control, may have profound biological functions. Positive feedback loops are ubiquitous in cellular regulation, with a major function to filter both the strength and temporal fluctuations of stimulating signals and to prevent inadvertent cell fate change. This network, however, may have an inherently slow response time, and the response is highly sensitive to the initial concentration of the substrate that lead to large cell-to-cell variation of temporal dynamics. This variation and slow dynamic may be problematic for processes such as neural crest formation and wound healing where precise and synchronized temporal control is crucial for generating collective responses of multiple cells. The expanded network shown in Fig. 5b allows transient increase of the initial substrate concentration, and solves the seemingly incompatible requirements for the simple positive feedback motif on robustness against fluctuations as well as fast and synchronized responses. It assures that despite a possible broad distribution of basal expression levels of the protein, cells are activated within a designated period of time at the presence of persistent activation signal, without sacrificing the filtering function of the feedback loop.
Cells use TOSS formed by a composite network to increase information transfer capacity Cells constantly encounter TGF-β signals with different strengths and duration, and must respond accordingly. It is well documented that biological networks reliably transmit information about the extracellular environment despite intrinsic and extrinsic noise in a subtle and functional way. However, quantitative analyses using information theory reveal that the dynamic of each individual readout is quite coarse with one or few bits. 9,10 This is a paradox. However, our results suggest that cells use multiple readouts to generate a TOSS with an expanded capacity to encode signal information and generate a far more subtle response system. For example, the SMAD motif has a refractory period due to the negative feedback loop and thus can accurately encode the duration information of TGF-β only within a limited temporal range. The GLI1 motif encodes information of longer TGF-β duration, which then saturates. This TOSS may be further expanded, such that the SNAIL1 motif itself possibly encodes information of longer TGF-β duration and relays to other transcription factors such as TWIST and ZEB, and leads to stepwise transition from the epithelial to the mesenchymal phenotype depending on the TGF-β duration. 34 Therefore although each motif has limited information coding capacity, a combination of motifs can code and transmit detailed signaling information. This is analogous to the design of a computer composed of many binary logic gates.
As with other signaling process, TGF-β signaling is context dependent, and the dynamic and regulatory network vary between cell types. 25,51 For the three cell lines we examined our results identify GLI1 as a major relaying factor for the TGF-β signaling. The inhibition experiments show that other possible peripheral links have minor contributions to SNAIL1 activation, while their weights may grow at time later than we examined. Consequently the present work has focused on the early event of TGF-β activation of SNAIL1, which is within 72 h for MCF10A cells. Nevertheless, the relay mechanism and the corresponding network structure identified here can be general for transmitting quantitative information of TGF-β and other signals. It is typical that an extracellular signal is transmitted through a canonical pathway with negative feedbacks and multiple non-canonical pathways, and these pathways crosstalk at multiple points, and Supplementary Fig. S6d gives some examples including IL-12, DNA double strand breaking, and LPS. Therefore, the mechanism revealed in this work is likely beyond TGF-β signaling.
Network temporal dynamics is a key for effective pharmaceutical intervention Upregulation of GLI1, and GSK3 and the responsive SMAD family has been reported in pathological tissues of fibrosis 52 and cancer, 49 and all three have been considered as potential drug targets. The present study emphasizes that in cell signaling timing is fundamental for function. The same network structure may generate drastically different dynamics with different parameters, as observed for different cell types. Consequently, effective biomedical intervention needs to take into account the network dynamics. We have already demonstrated that adding the inhibitors at different stages of TGF-β induction can be either effective vs. not effective on reducing SNAIL1 (by inhibiting pSMAD2/3), both effective (by inhibiting GLI1), and even opposite (by inhibiting GSK3). Actually, one may even exploit this dynamic specificity for precisely targeting certain group of cells while reducing undesired side effects on other cell types.
In summary through integrated quantitative measurements and mathematical modeling we provided a mechanistic explanation for how cells read TGF-β duration. Several uncovered specific mechanisms, such as expanding information transmission capacity through signal relaying, and reducing response times of positive feedback loops by increasing initial protein concentrations, may be general design principles for signal transduction.

TGF-β induce and inhibitor treatment
Cells for TGF-β induction and inhibitor treatment were seeded at~60-70% confluence without serum starvation. For TGF-β treatment, 4 ng/ml human recombinant TGF-β1 (Cell signaling) was added into culture medium directly. For inhibition experiment, 4 μM of LY2109761 (Selleckchem), 20 μM of GANT61 (Selleckchem), and 10 μM of SB216763 (Selleckchem) were used to inhibit SMAD, GLI, and GSK3, respectively. The medium was changed every day during treatment to keep the reagent concentration constantly. For reproducibility, we used cells within 10th-15th generations, same patches of reagents, serum, and tried to perform each group of experiments (e.g., those in Fig. 2c) together.

Immunofluorescence microscopy and data analysis
Cells were seeded on four-well glass-bottom petri dishes at~60% confluence overnight and treated with reagents (TGF-β1 and/or inhibitors). Three independent experiments were performed in every treatment. At designated time points, cells were harvested and stained with specific antibodies following procedure modified from the protocols at the Center of Biological Imaging (CBI) in the University of Pittsburgh. In general, cells were washed with DPBS for 5 min for three times followed by 4% formaldehyde fixation for 10 min at room temperature. Cells were then washed three times with PBS for 5 min every time. PBS with 0.1% TritonX-100 (PBS_Triton) was used for penetration. BSA of 2% in PBST was used for blocking before staining with antibodies. The first antibodies, anti-pSMAD2/3 (Santa Cruz), anit-SNAIL1 (Cell signaling), anti-GLI1 (Santa Cruz) were diluted by PBST with 1% BSA. Samples were incubated with the first antibodies at 4°C overnight. Then cells were washed three times with 10 min for each before being incubated with the secondary antibodies, anti-mouse Alexa Fluor 647 (Abcam), anti-rabbit Alexa Fluor 647 (Abcam), anti-gaot Alexa Fluor 647 (Abcam), anti-mouse Alexa Fluor 555 (Abcam), anti-rabbit Alexa Fluor 555 (Abcam), anti-gaot Alexa Fluor 555 (Abcam), anti-mouse Alexa Fluor 488 (Abcam), anti-rabbit Alexa Fluor 488 (Abcam), anti-gaot Alexa Fluor 488 (Abcam), for 1 h at room temperature. After antibody incubation, cells were washed with PBS_Triton for 5 min and stained with DAPI (Fisher) for 10 min at room temperature. Cells were washed three times with PBS_Triton for 5 min and stored in PBS for imaging.
Photos were taken with Nikon A1 confocal microscopy at CBI. The microscope was controlled by the build-in software, Nikon NIS Elementary. All photos, except the photo for GSK3 AA subcellular localization, were taken with plan fluor 40× DIC M/N2 oil objective with 0.75 numerical aperture and 0.72 mm working distance. The scan field were chosen randomly all over the glass-bottom area. For identifying the GSK3 AA subcellular localization, plan apo λ 100× oil objective with 1.45 NA and 0.13 mm WD was used. The 3D model of GSK3 AA overlapped with ERC and DAPI were reconstructed from 25 of z-stack images in 11.6 μm and videos were produced also by NIS Element software. To minimize photobleaching, an object field was firstly chosen by fast scan, then the photos were taken at 2014 × 2014 pixel or 4096 × 4096 pixel resolution, for generation large data or for photo presentation, respectively.
CellProfiller was used for cell segregation and initial imaging analysis as what described in Carpenter et al. 53 Image correction. To keep identical background through all images, background correction was performed before further image processing. For each image fluorescent intensities in space without cells were used as local background. Photos that have obviously uneven illumination and background fluorescence were removed from further processing. Otherwise the mean background fluorescent intensity was obtained through averaging over the whole image, and was deducted uniformly from the image.
Image segmentation. Cell number and position were determined by nuclear recognition with DAPI. The global strategy was used to identify the nuclear shape, and the Otsu algorithm was used for further calculation. Clumped objectives were identified by shape and divided by intensity. Next, using the shrank nuclear shape as seed, cell shape was identified by the Watershed algorithm. For identifying the clusters of GSK3 AA formed around a nucleus, the nuclear shape was shrank manually by 3 pixels and used as a new seed to grown the boundary with the watershed method until reaching background intensity level. All parameters were optimized through an iterative process of automatic segmentation and manual inspection.
Image quantification. Averaged fluorescence density and integrated fluorescence intensity were calculated automatically with CellProfiler. The amount of the GSK3 AA form was quantified as the sum of intensities of pixels belonging to the cluster formed around a nucleus. Concentrations of all other proteins were quantified by the average pixel intensity within the nucleus or cytosol region of a cell. Next, the quantified results were examined manually, and those cells with either cell area, nuclear area, or fluorescent intensity beyond five folds of the 95% confidence range of samples from a given treatment were discarded, which account for less than 1% of the cells analyzed. Immunofluorescence data were further processed and plots were generated using customized R codes and Matlab codes.
Quantitative PCR Cells were seeded in 12-well plastic bottom cell culture plates and treated as described above. Three parallel experiments were performed in every treatment. Total RNA was isolated with the TRIZOL RNA isolation kit (Fisher), and mRNA was reversely transcribed with the RNAscript II kit (ABI). The stem-loop method was used for microRNA reverse transcription. The qPCR system was prepared with the SYBR green qPCR kit with designed primers (Supplementary table 1) and performed on StepOnePlus real-time PCR (ABI).
Immunoprecipitation and silver staining Immunoprecipitation was performed with SureBeads magnetic beads (Bio-Rad) following a protocol modified from the one provided by the manufacture. We washed beads with PBS with 0.1% Tween 20 (PBS_Tween) for three times, then harvested cells by RIPA (Thermo) with proteinase and phosphatase inhibitor (Roche). Samples were pre-cleaned with 100 μl of suspended Protein G per 450 μl of lysis mixture. Antibodies Pathway crosstalk enables cells to interpret TGF-β J Zhang et al.
targeting GSK3 (Cell Signaling), GSK3 AA (Santa Cruz), and GSK3 D (Santa Cruz) were added into every 100 μl of bead mixture respectively. The mixture was rotated at 4°C for 3 h. Beads that were conjugated with antibodies were washed with PBS_Tween. An amount of 100 μl of precleaned lysis buffer was added into conjugated beads and rotated at 4°C overnight. Targeted proteins were eluded from beads by incubating with 40 μl 1× Laemmli buffer with SDS at 70°C for 10 min. For the samples an amount of 5 μl was used for western blot assay, and an amount of 30 μl was loaded for SDS-PAGE (Bio-Rad) and followed by silver staining (Fisher).

Network reconstruction and coarse-graining
The full network from TGF-β1 to SNAIL1 (Supplementary Fig. S1) was generated with IPA (Qiagen®). Specifically, all downstream regulators of TGF-β1 and upstream regulators of SNAIL1 in human, mice and rat were searched and added to the network. Then, direct or indirect relationships between every pair of regulators were searched and added to the network. After obtaining the whole network, regulators that have been reported to be activated later than SNAIL1 were removed. Examination of the network reveals that the network can be further organized into three groups: the TGF-β-SMAD-SNAIL canonical pathway, the TGF-β-GSK3-β-catenin pathway that has the most number of links, and others. We further noticed that GLI1 is a central connector of TGF-β, SMAD, GSK3 and SNAIL1. We performed western blot and IF studies on β-CATENIN and found that neither its concentration nor its location changes significantly on day 3, therefore we removed β-CATENIN from the network. In addition, previous studies report that the SMAD-GLI axis plays important role in TGF-β induced EMT. 31 Therefore we further grouped the network as the SMAD module, the GLI module, and the GSK3 module, as well as the remaining ones that we referred as "Others", and reached the network shown in Fig. 2a. Those molecular species not explicitly specified in Fig. 2a either have their effects implicitly included in the links, for example the link from TGF-β to GSK3, or are included in the links of "Others". This treatment is justified since our various inhibition experiments indeed showed that the three factors we identified affect SNAIL1 expression the most. These "other" species may contribute to snail1 activation at a time later than what considered in this work. Therefore we emphasize the network in Fig. 2a is valid only within the time window we examined, i.e., within 3 days after TGF-β1 treatment for MCF10A cells.

Data and code availability
The data and customized Matlab codes for simulation is available in following link: https://figshare.com/s/3d045c7da0db4fb4cd0f.