Cell signaling dynamics and transcriptional regulatory activities are variable within specific cell types responding to an identical stimulus. In addition to studying the network interactions, there is much interest in utilizing single cell scale data to elucidate the non-random aspects of the variability involved in cellular decision making. Previous studies have considered the information transfer between the signaling and transcriptional domains based on an instantaneous relationship between the molecular activities. These studies predict a limited binary on/off encoding mechanism which underestimates the complexity of biological information processing, and hence the utility of single cell resolution data. Here we pursue a novel strategy that reformulates the information transfer problem as involving dynamic features of signaling rather than molecular abundances. We pursue a computational approach to test if and how the transcriptional regulatory activity patterns can be informative of the temporal history of signaling. Our analysis reveals (1) the dynamic features of signaling that significantly alter transcriptional regulatory patterns (encoding), and (2) the temporal history of signaling that can be inferred from single cell scale snapshots of transcriptional activity (decoding). Immediate early gene expression patterns were informative of signaling peak retention kinetics, whereas transcription factor activity patterns were informative of activation and deactivation kinetics of signaling. Moreover, the information processing aspects varied across the network, with each component encoding a selective subset of the dynamic signaling features. We developed novel sensitivity and information transfer maps to unravel the dynamic multiplexing of signaling features at each of these network components. Unsupervised clustering of the maps revealed two groups that aligned with network motifs distinguished by transcriptional feedforward vs feedback interactions. Our new computational methodology impacts the single cell scale experiments by identifying downstream snapshot measures required for inferring specific dynamical features of upstream signals involved in the regulation of cellular responses.
Cells continuously sense a variety of physical and chemical signals and respond suitably to changes in their environment [1–3]. Depending on the nature of the environmental change, a subset of membrane receptors gets stimulated, eliciting specific signaling pathways and enabling cells to make informed decision downstream [4]. Recent studies indicate that the functional information to trigger a specific downstream response is typically carried through activity change in one or more signaling molecules, transcription factor activities or protein substrates [5–7]. Considering the multifunctionality of signaling pathways [8]), it has been argued that this mode of information transfer is more rich and complex than encoding it through merely binary on/off steady states [9–12]. For example, depending on the cell type, transient activation pattern of stimulated extracellular signal-regulated kinases (ERK) leads to cell proliferation, whereas sustained pattern of ERK may lead to cell differentiation [13]. However, the mechanism through which the functional information encoded within the temporal activity change is processed and decoded to very fine alterations in an upstream signaling pattern still remains unclear. Instantaneous activity change between upstream and downstream events has been typically considered to be the natural mechanism of information processing, thus, limiting our understanding of dynamical patterns [14–17]. Therefore, a challenging but relevant task is to elucidate the mechanism through which features or properties of signaling dynamics (time dependent activity-change in signals) encode information, and how cells decode this information via complex signaling pathways with high specificity, when variability is the norm. Mounting evidence from single cell microscopy [18], flow cytometry [19], single cell PCR [20], single cell RNA-seq [21,22], and single cell mass cytometry [23] studies reveals that variability in post-transcriptional or post-translational events, even within homogeneous condition, is ubiquitous [24–26]. Such results has generated considerable interest in employing single cell resolution data to elucidate the role of variability and characterize underlying mechanism of intracellular information processing [3,27]. Since comprehensive network motif identification studies have unveiled a strong interdependence between dynamics, motif structure, and specific function [28], role of network topologies during information processing also needs to be considered [29]. For instance, feedforward loops have an ability to generate a transient activity or protect against brief fluctuations depending on the nature of their interactions [28], whereas positive feedback loops induces a response delay and hysteresis in the network dynamics [30]. Negative feedback mechanisms protect the information transfer by saturating the signaling pathways from upstream fluctuations [31]. In sum, precise orchestration of cellular events is a consequence of an efficient signaling modulation code involving upstream signaling dynamics and network motifs [1,32,33]. Hence, we hypothesize that functional information to direct cells towards specific responses is encoded in the temporal features of the signaling dynamics, and that the regulatory network motifs facilitate decoding this information. We test our hypothesis through analysis of an integrated model of signaling dynamics and gene regulatory networks. Both gene regulatory networks and communication networks have a similar underlying operating principle [34,35]. As shown in Fig 1 the goal is to encode relevant information in the characteristic features of the input signals and transmit them over a channel so they can be decoded at the receiving end [17,36,37]. Molecular signaling pathways elicited by receptor activation in cells which are involved in the regulation of physiological responses such as growth, differentiation and programmed apoptosis, among others [11,38], have been shown to exhibit such properties [3,9]. Unlike communication networks, there are limitations to the predictability of downstream responses in a biochemical network for any given input signal, and in particular at single cell levels [39]. This is because it is difficult to infer the strength of an input stimulus received by a single cell based on the instantaneous downstream responses alone. The goal of this manuscript is to assess how variation in dynamic signaling pattern influences downstream regulation at the single cell level. Identifying the mechanism of how a specific downstream phenotype could be triggered by different combinations of kinetic modulation in an upstream signaling pattern is critical to our understanding of the cellular decision making process. Here we employ a mechanistic ODE based model [40], that is comprised of multiple receptor stimulated signaling kinase cascades activating immediate early response genes (IEGs) (c-Fos and c-Jun), and transcription factors (TF) (AP-1), which turn on late response target genes (right side of Fig 1) [41]. This two-tier regulatory network is comprised of a feedback as well as a feedforward transcriptional network motif, making it a suitable candidate to test our hypothesis. We characterized a set of independent dynamical features for each input signal (ERK, FRK and JNK in Fig 1) through a novel phenomenological model (see Fig 2A for graphical illustration). Variance based global sensitivity analysis was used to study how uncertainty in signaling features affects uncertainty in downstream IEG and TF responses, and to also quantify the relative importance of each feature at every timepoint [40,42]. Information theory was used to measure mutual information between signaling feature and the downstream gene regulatory responses. Mutual information quantifies the amount of information that is carried by a signaling feature, and allow us to evaluate its capacity to make a binary transcriptional decision downstream. We further use decision tree analysis [43], a predictive modeling technique, to identify decision rules for signaling features that significantly shape the IEG and TF response patterns. As single cell studies become more advanced, it has become more evident that cells may use dynamics of signaling species to encode important information that can be decoded to generate distinct cellular responses. We use our earlier published single cell gene expression data [20] for specific transcriptional target outputs to identify the upstream signaling features that are most likely to be transmitted into specific phenotypes via generation of distinct target gene expression levels. Our results revealed that while signaling dynamics have many potential functional information encoding features, transcriptional regulatory network motifs play a significant role in limiting which features are significant at different phases of the stimulation. Moreover, estimating the information transferred by each signaling feature showed that in the majority of cases a single feature is insufficient to lead to a binary transcriptional decision (i. e. individual features carried less than one bit of information). In this light, decision tree analysis demonstrated that a combination of signaling features constitutes a feature based modulation code involving crosstalk in order to drive the cellular decision making. Given extensive single cell variability, our methodology is of practical utility to experimentalists who use single cell resolution data and are typically confronted with the variable aspects of the cellular activities [44,45]. By focusing on dynamical features of the inputs to a signaling system, we have provided a novel framework which is broadly applicable to many signaling systems.
A phenomenological model for dynamics of a transient signaling kinase was developed to characterize independent time-invariant signaling features based on dynamical phases of delay, activation, peak retention and deactivation kinetics (Fig 2A; Methods section). An exponential ascent or decay model for temporal phases was derived to establish duration and time-constant features (see Table 1 for details). In order to capture observed variability in signaling dynamics, we introduced simultaneous variation in all of the identified upstream signaling features using a differential variation schema (Table 1 and S1B Fig). Our results revealed two distinct aspects of signaling variability. First, we observed that simultaneous variation in signaling features captured the variation of the entire signaling profile (shown as gray region in the Fig 2B, 2C and 2D). Second, these variations were found to be within experimental range measured from bulk tissue samples, shown as red error bars in the plots. Our results suggest that the variability in signaling features spans the combined effects of intrinsic and extrinsic fluctuations. A combination of three distinct dynamical profiles for each signaling kinase ERK, FRK and JNK, respectively, was selected at random and simulations were performed for each of these combinations. A total of 100,000 simulations were conducted. The response of immediate early genes (IEGs) and transcription factors (TFs) resulting from these simulations revealed that variability in dynamical features propagates downstream as signal progresses in the network (gray regions in Fig 2E, 2F, 2G, 2H and 2I). Specifically, variability in downstream AP-1 TF, the aggregate of the two dimers, matched well with the observations from cell population measures [46], shown in red error bars, implying that variations in signaling features captures the biological variability at different levels in the network. Each dynamic signaling feature impacted the overall downstream response both positively and negatively, depending on the nature of the variations and dynamical interactions involved. Downstream response to receptor stimulation was divided into three phases: (1) Early phase where cells predominantly utilized the steady state TF proteins as well as triggered new IEG expression, (2) Intermediate phase where cells utilized new proteins transferred from IEG expression to modulate AP-1 TF composition and levels, and when the heterodimer ppc-FOS: ppc-JUN was the dominant contributor to AP-1 aggregate [47,48], and (3) Late phase where ppc-JUN: ppc-JUN homodimer was the primary contributor to AP-1 TF [49]. It is important to note that transcription of c-Jun gene is slower in comparison to that of c-Fos gene [40], however, the translation process of c-Jun mRNA is much faster in comparison to that of c-Fos mRNA [41]. Variability in downstream TF dynamics (shown as density plots in Fig 3A) were decomposed to quantify the contribution of each signaling feature at different phases of the stimulation. Variance based global sensitivity analysis was performed to measure these contributions as first order and total sensitivities (see Methods). These evaluations revealed that AP-1 dynamics (both aggregate and individual dimers) were sensitive to only a select subset of signaling features at any given time (Fig 3B). Transcription regulatory pathways contain a small set of recurring regulation patterns, called network motifs that are critical in processing of signaling information and in shaping of the response pattern. Our network contains two such recurring network motifs (S2 Fig): (1) FRK/ERK module showing a coherent feed forward interactions between receptor stimulated ERK pathway to activate ELK1 and subsequently transcribe IEG c-Fos, which is further activated (phosphorylated) by FRK to produce the heterodimer. (2) JNK module represents a positive feedback loop by expressing IEG c-Jun buffered via pp-ATF2: ppc-JUN dimer complex, and activation (phosphorylation) of the homodimer. Estimating the lower bound of the Mutual Information (MI) between a signaling feature and downstream TF response revealed that majority of the information carried by individual features to TFs were less than one bit (See Methods and Fig 4). In this context, MI corresponds to the decoding of the TF level which can discriminate between the two levels of a dynamic signaling feature. Analysis of the information transfer at 20 minutes (early phase) revealed that α2F and τ1F transferred nearly equal amount of information (∼ 0. 4 bits) to both heterodimer as well as to AP-1. In contrast, at 40 minutes (intermediate phase) and 60 minutes (late phase) the amount of information transferred by activation time-constant of JNK (τ1J) to AP-1 was the highest (∼ 0. 4 and ∼ 0. 45 bits, respectively). However, τ1J carried ∼ 0. 7 and ∼ 0. 5 bits of information to homodimer at 40 and 60 minutes, respectively, a complete reversal to the pattern observed in the Total AP-1. Importantly, the time-constant feature of FRK (τ3F) carried ∼ 0. 4 bits at 40 minutes and ∼ 1. 2 bits at 60 minutes, suggesting that the network has the capacity to decode this feature to make binary decisions at the later phase of the dynamic response. In contrast, the FRK time-constant τ3F carried less than 0. 1 bit of information to Total AP-1. The signaling features of ERK had minimal information transfer to the TFs, in consistent with the network topology, as the effect of changes in ERK on the TFs is indirect and conditional on the dynamic features of FRK and JNK signaling. We extended our results on global sensitivity analysis and information transfer to study whether the variable responses at the single cell level could be analyzed for inferring the dynamic features of upstream signaling. We considered a typical case involving snapshot measures from single cells, e. g. , gene expression levels, or immunocytochemistry data on activated TF protein levels [51]. Decision tree analysis was performed (see Methods) on downstream TF activity phenotypes at three different timepoints: 20 minutes (S5 Fig), 40 minutes (S6 Fig) and 60 minutes (Fig 5). We observed that the TF activity response at any given time point was a distribution, and a function of varying upstream signaling features. Three distinct phenotypes at each time point were defined as High, Mid, and Low TF levels (Fig 5). The High and Low phenotypes represent the upper and lower margins of the distribution, respectively, and the Mid phenotype corresponds to the median of the distribution. We considered the TF phenotypes of Total AP-1 (Fig 5A), heterodimer ppc-FOS: ppc-JUN (Fig 5E), and homodimer ppc-JUN: ppc-JUN (Fig 5I). Mapping these TF phenotypes to corresponding kinase profiles resulted in distinguishable signaling features that would yield distinct TF response dynamics (Fig 5B, 5F and 5J). The combination of key features varied depending on the specific TF phenotype analyzed: JNK features appear to be distinct for the Total AP-1 and homodimer phenotypes (Fig 5B and 5J), whereas FRK features were specific to the stratification of the heterodimer phenotypes (Fig 5F). We developed decision trees to organize these signaling features into a hierarchy of rules defining the ranges of variation that stratify activation profiles of downstream TF phenotypes. For the Total AP-1 TF snapshot response measure at late phase, the decision tree analysis revealed a hierarchy of rules based on three signaling features of JNK (τ1J, τ2J and τ3J) as primarily yielding downstream phenotypes (Fig 5C). Our results revealed that fast activation of JNK (τ1J < 5. 3) and very slow peak deactivation of JNK (τ2J > = 25) were highly likely to drive cells to a High AP-1 phenotype, even in presence of variations in other kinase signaling features. In contrast, slow activation of JNK (τ1J >= 5. 3), fast peak-deactivation of JNK (τ2J < 16), and fast deactivation of JNK (τ3J < 22) were very likely to drive cells to Low AP-1 phenotype. We assessed the decision tree analysis results in an alternative visualization layout considering the top two dominant features from the decision trees in a scatter plot, and mapped the rules as well as the phenotype categories to the plot to observe how these features separated cells with a specific phenotype (Fig 5D). The extremes, High and Low AP-1, were well separated by specific intervals of τ1J and τ2J, with faster activation and slower inactivation leading to a High AP-1 phenotype. In the case of heterodimer phenotypes, τ3F was found to be the only dominant feature separating all of the three phenotypes (Fig 5G and 5H). On the contrary, the homodimer phenotypes were segregated by a combination of τ1J, τ2J, τ3J and α1J (Fig 5K and 5L). The results from decision trees were in agreement with those from analysis of MI. Features which carried one or more than one bit of information to downstream TF responses were sufficient to individually distinguish High, Mid, and Low phenotypes. For example, α1J separated homodimer phenotypes at 20 minutes (S5K and S5L Fig), and τ3F separated heterodimer phenotypes at 60 minutes (Fig 5G and 5H). α1J carried ∼ 0. 85 bits of information to homodimer at 20 minutes, and τ3F carried ∼ 1. 2 bits of information to heterodimer at 60 minutes (Fig 4). In summary, decision tree analysis revealed that an individual feature which transferred ≥ 1 bit of information corresponds to the distinct downstream TF phenotypes, with potential consequences for fine tuning cellular decision making in response to external stimuli. We further analyzed bivariate gene expression patterns of immediate early genes (c-Fos and c-Jun) as observed in single cell experiments (Fig 7A, data is from [20]). Despite the cells being picked from a pool of neurons of same type under homogeneous conditions, experiments show that the variability in the gene expression levels is extensive [53] (almost 4-8 fold variation in the abundance). In order to identify the key upstream signaling features that lead to distinct phenotypes of the observed response distributions (Fig 7A), we employed a model based approach to analyze statistically similar distribution of the single cells from our simulation results. Four groups were identified based on the high and low 25th percentile of the marginal distributions: (1) high c-Fos, high c-Jun, (2) high c-Fos, low c-Jun, (3) low c-Fos, high c-Jun, and (4) low c-Fos, low c-Jun. We found that the key controlling features established from our earlier analysis (Fig 6), separated the four IEG groups distinctively in the signaling feature space (Fig 7B). Decision tree analysis (Fig 7C) revealed that: (1) slow rate of ERK peak deactivation (high τ2E) and slow rate of JNK peak deactivation (high τ2J) produced high c-Fos and high c-Jun gene expression pattern, (2) slow rate of ERK peak deactivation (high τ2E) but fast rate of JNK peak deactivation (low τ2J) produced high c-Fos and low c-Jun gene expression pattern, (3) fast rate of ERK peak deactivation (low τ2E) but slow rate of JNK peak deactivation (high τ2J) produced low c-Fos and high c-Jun gene expression pattern, and (4) fast rate of ERK peak deactivation (low τ2E) and fast rate of JNK peak deactivation (low τ2J) produced low c-Fos and low c-Jun gene expression pattern. Our results indicate that the variability is the gene expression patterns in single cells is likely to contain information on features of the upstream signaling dynamics. Our analysis thus far revealed that specific features of signaling dynamics correspond to the patterns of downstream gene regulation and TF activity. The results also indicated that distinct signaling features can be decoded depending on the downstream location in the regulatory network. This motivated us to consider whether there was a higher order organization of which features being decoded and where in the network. We performed exhaustive global sensitivity analysis and information transfer assessment to consider each node in the network as the output of interest. For each node, we computed the sensitivity indices corresponding to the dynamic signaling features as well as the information transfer between the signaling feature at the corresponding node levels. We represented the results as a Sensitivity Map and Information Transfer Map. The Sensitivity Map lays out all of the controlling features (Si > 0. 1) for different nodes in the regulatory network and at different phases of the stimulation (early, intermediate and late) (Fig 8A). Similarly, the Information Transfer Map lays out the extent of information transferred by each feature to different network decoders (mutual information between signaling features and the downstream responses) and at different phases of the stimulation (Fig 8B). Unsupervised clustering of these maps yielded two major groups corresponding to the two regulatory motifs that drove distinct downstream activity patterns: transcriptional feedforward and transcriptional positive feedback interactions (Fig 8C). Both of these maps independently separated regulatory motifs along the network decoders thereby revealing the interdependence of information processing by signaling features and by network topology.
We present a new paradigm in which dynamic signaling features and not molecular abundances are introduced as the carriers of functional information. We developed a novel empirical model for signaling profiles based on transient phases of activation, peak retention and deactivation kinetics to characterize the dynamical aspects of the signal. We then performed dense Sobol sampling over an uniform high-dimensional space of signaling features to map the functional space of signaling dynamics which spans all possible essential dynamical profiles for three upstream signaling kinases, namely ERK, FRK and JNK respectively. We employed variance based global sensitivity analysis, estimated mutual information, and developed decision trees to identify the key controlling features. The amount of the information carried to various network components were also analyzed, specifically to the IEG and TF responses. Our results showed that the key controlling feature must carry at least one bit of information in order to distinctively separate downstream IEG and TF dynamical response patterns. Our analysis revealed that in the majority of cases an individual signaling feature carried insufficient amount of information (less than one bit) to directly lead a binary transcriptional regulation, and therefore a combination of features influenced downstream regulatory process. Cellular decision, hence, is a consequence of decoding of multiple dynamic signaling features, rather than a direct function of the abundances of functionally active signaling effectors. With recent advancements in single cell technologies such as fluorescent labeling techniques, optogenetics, gene-expression profiling, and time-lapse microscopy, it is evident that cells transduce functional information through time-dependent activity changes in one or more signaling molecules [9–11,54–56]. For example in a recent single cell quantitative live-cell imaging study, it was found that individual macrophages exhibited analog responses to LPS-induced NF-κB dynamics [57], in contrast to the digital responses previously reported in cytokine treated non-immune cells [58,59]. Cells, therefore, are capable of interpreting a stimulus in a continuous manner for which the time-dependent aspect of the signal represents an embedded functionality. We explored a scheme where the cell’s signaling code involves dynamic signaling features. These features represent different kinetic modulation of the signal as the functional aspect of a molecular activity change. Simultaneous variation in these features yielded signaling dynamics as variable as observed experimentally, and also bounded within the limits of the population measure. As the signal progressed in the network, the variability through dynamical features also progressed resulting in a non-uniform response distribution downstream, in consistent with other studies [24,26,60,61]. Our results showed that this downstream response distribution were also within the biological variability reported in the population response. Using global sensitivity analysis, decomposition of the variance of the downstream IEG and TF responses revealed that the response distributions were essentially a result of multiplexing of various signaling features. Our results implied that pleiotropic regulators such as signaling pathways are likely to employ a feature-based modulation code to transduce intracellular information and elicit context-specific responses. This code includes combinations of time-invariant signaling features including duration, time constant, amplitude and/or frequency of dynamics to ensure specificity in the downstream responses. Further, our approach employs decision tree analysis to identify the rules involving upstream signaling features to understand the nature of decision making, enabling us to unravel the complexity of cellular information processing. Our approach points to a novel modulation code involving dynamic features of signaling, through which the gene regulatory phenotypes may be shaped. In such a formulation, the manipulation of signaling towards desired gene regulatory responses involve alterations in the signaling onset, duration, etc [62]. This is in stark contrast to an abundance based code that would use over/under-expression or knock in/out based interventions. Analyzing the implications of transcriptional variability at the single cell level must include the interplay between the properties or features of the upstream signaling dynamics [41,60]. Decision tree analysis allow us to establish the functional implications of variability in signaling dynamics in cellular decision and elucidate how this variability may result in cellular (or phenotypic) heterogeneity. Our results revealed that correlations between activity levels of signaling molecules and downstream responses are insufficient to establish a physiologically relevant code. For instance AP-1 activity levels were found to be sensitive to a different subsets of signaling features at different phases of the stimulation. Our analysis found that the regulatory network is as essential an element in cellular decision making. Intracellular information decoding is equally as critical as information encoding. Similar observations have been reported in a single cell study where dual reporters were used to characterize live-cell response through LPS induced NF-κB dynamics, and to learn how they correlated with inflammatory cytokine gene-expression output [57]. When the instantaneous ratio of abundance levels was calculated, the transcriptional regulator (NF-κB) and downstream transcription (TNF-α gene reporter) were poorly correlated. An unbiased regression analysis of NF-κB dynamics over entire stimulation time course was found to be the key determinant of the transcriptional output, suggesting the limitation of comparing molecular abundances. In the same study, it was observed that the macrophages showed a robust discrimination to different LPS doses. This was because the information decoding involved a feedback dominance switching which stratified the immune response despite the limited information encoding capacity of the signaling network [17]. Hence, our dynamic feature based information processing framework can be readily employed to study cellular decision making in a broad range of biological contexts. Variability in signaling dynamics [63] and in transcriptional regulation [25,41] have been typically studied in isolation of each other. For instance, studies on the dynamics of ERK2 translocation in single human cells reveal a large variation in basal ERK2 nuclear levels [63,64]. Besides peak fold change, other features of the dynamic response pattern such as final fold change, the peak delay and timing of activation also exhibit high variability [65]. Similar variability is also seen among the single cell time courses of NF-κB dynamics treated across the range of LPS concentrations [57]. On the other hand, the seminal studies of transcriptional variability in bacteria and yeast have demonstrated that inherent stochastic noise in gene expression (intrinsic variability), as well as variability in cellular transcription factor activity levels (extrinsic variability), can result in cellular heterogeneity [26]. We have also reported similar observations in mammalian cells where cell-to-cell variability of transcriptional change was shaped by the inputs to the cells, suggesting the patterns of response distribution are based on a systematic processing of inputs instead of being a noise around the mean [20]. Our framework bridges the variability between the signaling dynamics and transcriptional regulatory responses through a model based approach. Most interestingly, our results demonstrate that even a single snapshot response measure of single cell transcriptional activity can be informative of upstream signaling dynamics. Further analysis revealed that different network components (IEG vs TF) decode different aspects of the signaling dynamics, indicating that a network component can only decode a portion of the signaling history. For instance, AP-1 response dynamics were informative of a combination of activation and peak retention kinetics of the JNK signal, whereas c-Fos and c-Jun IEGs were informative of peak deactivation kinetics of ERK and JNK signals, respectively. This distinction in dynamic decoders is a result of multiplexing of signaling features via regulatory network motifs. Feedforward interactions imposed a bottleneck to information processing, where early phase features such as delay and activation kinetics of FRK did not transmit information to downstream TFs. In contrast, positive feedback loop introduced a time lag in information processing where the late phase of the downstream AP-1 TF was informative of activation kinetics (early phase) of the JNK signal. Other studies on signaling dynamics, particularly those considering networks with feedback and feedforward loops, use much longer term analysis than just one hour [66,67]. Opportunities exist for longer-term dynamics to unravel additional complexities of these regulatory network motifs which may not be captured by shorter-term dynamics. Studies to date have analyzed the mechanism of cellular decision making as related to dynamics of signaling activities by evaluating the correlation between the upstream and downstream molecular abundances. For instance, high correlation between EGF binding to EGFR and subsequent cell proliferation [64], or high correlation between increase in IL2 and IL4 receptor trafficking and consequent T cell proliferation [68,69] have been interpreted as information transfer from the activated receptor levels to the downstream effector activities. However, inferences based on instantaneous correlations do not consider how cells process the time-dependent activity changes. The implicit assumption to these studies is that instantaneous change in molecular abundances between the upstream and downstream events is the natural mechanism of cellular information encoding [31,70]. Such an approach may establish a causality between the two domains, but limits our understanding of any functional capabilities that could be gained by processing dynamical pattern. Our approach of reformulating the problem overcomes these limitations for a novel analysis of cellular information processing that takes dynamics into account. Emerging single cell studies are illustrating the use of high temporal resolution data sets to unravel the mechanism underlying cellular decision making process. For instance in a recent study, a microfluidic setup was used to dynamically administer inhibitor treatment to alter the amplitude, frequency, and duration of nuclear localization of transcription factor Msn2, while a fluorescent reporter of Msn2 transcriptional activity was used to detect the downstream response [7,33]. The study revealed that different expression patterns correlated with distinct dynamical features of administrators. As another example, a short duration of LPS stimulus does not elicit an immune response, whereas persistent LPS stimulus yields an innate immune response [71]. Such results are consistent with our findings that the information transfer underlying cellular decision making involves dynamic aspects of signaling. We found that the dynamic decoding is in close alignment with the network topology, as revealed by unsupervised clustering of the newly developed Sensitivity and Information Transfer Maps. Our findings suggest a novel dynamic modulation code, where signaling features in coordination with the network motifs shape the downstream regulatory patterns. Our approach provides a novel framework to study the structural features, i. e. , static aspects of a network, in integration with signaling and transcriptional regulatory kinetics, i. e. , dynamic aspects of the function, to understand information transfer underlying cellular decision making process. Understanding these dynamic encoding-decoding principles would enable sophisticated strategies to manipulate signaling dynamics to better control cellular responses in a broad range of biological contexts.
We defined gene regulatory model as a combination of regulatory network ℱ and external input signal fs as y ˙ = F (t, y 0, k, f s (X, t) ) (1) where species y ∈ ℝ s, parameters k ∈ ℝ m and signaling features X. Activation of receptor due to ligand binding on the cell surface elicits and localizes multiple signaling kinases [72]. In our model, extracellular-signal-regulated kinase (ERK), c-FOS regulating kinase (FRK) and c-Jun N-terminal kinase (JNK) were activated in cells through receptor stimulation, to transcriptionally and post-translationally regulate the activator protein-1 (AP-1) family of TFs (See [40,73] for details on model development). Family of AP-1 TFs considered here are: the heterodimer ppc-FOS: ppc-JUN and the homodimer ppc-JUN: ppc-JUN. These activated AP-1 TFs subsequently bind to promoters of key target genes resulting in gene expression changes [20,74]. Activation of target genes is involved in regulation of blood pressure and the development of hypertension in brainstem [74–76]. Hence, AT1R stimulated multiple signaling kinases in brainstem neurons are tightly coupled with the specific physiological function of controlling blood pressure. This mechanistic model was implemented as a set of ordinary differential equations that integrated three signaling kinases viz. FRK, ERK and JNK with activation of AP-1 TFs (Fig 1). Experimentally measured time course activity change in three kinases were interpreted as input signals [40,46], and subsequent changes in the activity of IEGs and TFs as the downstream responses [40,77]. Details of the model equations are reported in S1 and S2 Tables. The model code in MATLAB format is available through accession No. 185122 on the ModelDB resource [78] at https: //senselab. med. yale. edu/ModelDB/showModel. cshtml? model=185122. In our model, signaling kinase profiles were assumed to be an outcome of enacting activation, peak retention, and deactivation processes, initiated in parallel by receptor stimulation through ligand binding (input stimulus) [79,80]. A simplified phenomenological interpretation of the phosphorylated kinase signaling profile can then be developed as the difference between first order activation and peak processes, and both of them opposed by a first order deactivation processes, represented as: f s (t) = (z 1 (t) - z 2 (t) ) - z 3 (t) (2) where fs is the dynamics of the signaling kinase, z1 are the activation processes, z2 are the peak processes, and z3 are the deactivation processes initiated by a receptor stimulus as shown in control diagram (S1A Fig). z 1 ˙ (t) = k a z 1 u (t) - k d z 1 z 1 (t) (3) z 2 ˙ (t) = k a z 2 u (t) - k d z 2 z 2 (t) (4) z 3 ˙ (t) = k a z 3 u (t) - k d z 3 z 3 (t) (5) where u (t) is a function representing the upstream stimulus, k a z i is the activation parameter of the dynamical process zi and k d z i is the deactivation parameter of zi. To characterize the dynamical signaling features of the activated kinase signal, a transfer function model can be derived by taking the Laplace transformation on each of these differential equations. Considering all three processes as analogous to each other in their behavior yields the following form. L { z ˙ (t) } = L { k a z u (t) } - L { k d z z (t) } (6) s z ^ (s) = k a z u ^ (s) - k d z z ^ (s) (7) F (s) = z ^ (s) = (k a z / k d z (1 / k d z) s + 1) u ^ (s) = (λ τ s + 1) u ^ (s) (8) For u (t) as a heavyside unit step function, the inverse laplace transformation (ℒ − 1 { F (s) }) results in an exponential ascent or decay model depending on the process [81]. Note that each of the activation, peak, and deactivation dynamical processes are defined with specified delays, αi, with a following constraint; α1 + α2 + α3 ≤ (tf − ti), where α1 is the delay in activation processes, α1 + α2 is the delay in peak processes, α1 + α2 + α3 is the delay in deactivation processes, ti is the initial time point and tf is the final time point. Resulting dynamics of the signaling kinase typically follow a transient signaling pattern over minute timescales [79,82]. This transient signaling profile, represented by fs in Eq (9), was observed to have four phases of dynamics: delay, activation, peak retention and deactivation phases. The mathematical form for last three phases are represented as a function of two independent features: (1) a duration feature, α and (2) a time constant feature, τ, as shown in Eq (9). The resulting function of signaling dynamics is continuous and differentiable within a given boundary condition (from ti = 0 to tf = 60 minutes in our case) (See Fig 2A for a graphical representation). f s (X, t) = λ { 0, t ≤ α 1 Delay 1 - exp (- (t - α 1) τ 1), α 1 < t ≤ (α 1 + α 2) Activation exp (- (t - α 1 - α 2) τ 2), (α 1 + α 2) < t ≤ (α 1 + α 2 + α 3) Peak Retention exp (- (t - α 1 - α 2 - α 3) τ 3), (α 1 + α 2 + α 3) < t ≤ 60 Deactivation (9) where X ∈ {α1, τ1, α2, τ2, α3, τ3} represents the collective set of all signaling features of a kinase signal. The signaling dynamics of each activated kinase is represented by six independent features, making it a total of eighteen for three kinases ERK, FRK and JNK, respectively. List of all independent features, their least square fitted values with their appropriate experimental observations, and the range of induced variation to each signaling feature are reported in Table 1. To computationally depict the variability of single cells, duration features which are in time domain were varied up to a multiplicative factor of 2 to 4, whereas time constant features were varied up to two-fold. In this study, a biochemical regulatory network is considered as the channel which processes various dynamical aspects of the signaling pattern and generates a downstream response pattern (Fig 1). Because multiple interactions occur simultaneously in a signaling pathway, a true mathematical relationship between the dynamical signaling features and a downstream response is difficult to construct. Hence, to evaluate the influence of signaling features on a downstream response, we turn to global sensitivity analysis. Global sensitivity analysis investigates how uncertainty in key signaling features of a biochemical model affects uncertainty in its response and quantify the relative importance of these features. We here use Sobol’s high dimensional model representation (HDMR) of biochemical regulatory networks to perform variance based global sensitivity analysis and to quantify the influence of signaling features [42,83]. The gene regulatory model is considered as a multivariate function of signaling features X ∈ {α1, τ1, …, αn, τn}, where these features acts as inputs, influencing the output IEG or TF responses, ℱ, either in an independent and/or cooperative way as shown below. F (X) = f 0 + ∑ i ≤ 2 n f i (X i) + ∑ 1 ≤ i < j ≤ 2 n f i j (X i, X j) + ⋯ + f 1,2, ⋯, 2 n (X i, X 2, ⋯, X 2 n) (10) where ℱ (X) ≡ ℱ (t, y 0, k, f s (X, t) ) as shown in Eq (1), f0 denotes the mean value of the ℱ (X) over the entire variation range of X. The first order function fi (Xi) describes the independent behavior of Xi on output ℱ (X). In addition, the functions also preserve a nonlinearity which is essential for establishing a general case. The second order function fij (Xi, Xj) represents the cooperative behavior between two input features Xi and Xj on output ℱ (X). Similarly, 2nth order functions represents the cooperative effects of 2n signaling features acting together to influence the output ℱ (X). We began this analysis by considering variability in signaling dynamics as a function of variability in signaling features, αi and τi in Eq (9), as a result of continuous intrinsic and extrinsic fluctuations in cells [25]. The simultaneous Monte-Carlo sampling for each signaling feature was employed to induce variability, in consistent with experimental observation, using function ‘sobolset’ and ‘scramble’ in the Matlab computational environment [84] (Table 1 for details on range of variation). HDMR decomposition measures the individual as well as higher order cooperative measurements of the variations in output (downstream response) in terms of variations in input (s) (or signaling feature (s) ). F (X) ≈ ∑ i ≤ 2 n f i (X i) ︸ first order + (∑ 1 ≤ i < j ≤ 2 n f i j (X i, X j) + ⋯ + f 1,2, ⋯, 2 n (X i, X 2, ⋯, X 2 n) ) ︸ higher order (11) Total variance V[R] for a response output R ≡ ℱ (X) can be decomposed into individual variances of each HDMR component function using mutual orthogonality [83,85]. The response R is the activity level of the downstream IEG or TF under observation. This is a very useful property of HDMR in determining how the uncertainty in downstream response (R) is influenced by the uncertainties in signaling features (X). V [ R ] = ∑ i ≤ 2 n V (E [ R | X i ]) + ∑ 1 ≤ i < j ≤ 2 n V (E [ R | X i, X j ]) + ⋯ + V (E [ R | X 1, X 2, ⋯, X 2 n ]) (12) When variation in each feature is assumed to be unbiased (i. e. Xi’s having a uniform distribution) S1 Fig, the Sobol sensitivity indices are defined as the ratio of decomposed partial variances to the total variance as shown below, where Si represents first order sensitivity indices, and S i T represents total sensitivity indices. S i = V (E [ R | X i ]) V [ R ] (13) S i T = V (E [ R | X i ]) + ∑ i j V (E [ R | X i, X j ]) + ⋯ + V (E [ R | X 1, X 2, ⋯, X 2 n ]) V [ R ] (14) The difference between the two S i T − S i is considered as the higher order sensitivity index. The conditional variance in the sensitivity indices were estimated using Saltelli’s approach [86]. While estimating the indices, we assume that the biochemical system responded sufficiently and quickly to reach steady state before there was any change in the signaling feature (otherwise the system would still be responding to a previous fluctuation in the features when the next change happened). The objective of the global sensitivity analysis is to identify the most important signaling feature affecting a downstream response at a given time. This is determined through identification of a feature which, when varied simultaneously along with all the features, leads to the greatest reduction in the variance of the downstream response. Similarly, second most important feature can be determined which leads to second greatest reduction in the variance of the downstream response, and so on until all features are ranked in the order of importance. Estimation of sensitivity indices captures this reduction in the variance through a continuous metric, bounded between 0 and 1 [86]. A sensitivity index of zero means that the associated feature is non-influential to a downstream response, while the higher the index the more influential the feature. If the sensitivity indices were lower than a preestablished threshold value (in our case Si < 0. 1), we assumed that the signaling feature was not sensitive enough to lead to significant reduction in the response variation and was not important. This criteria was used throughout the study and in generation of the heatmap in Fig 8. We applied information theory to quantify the amount of information transduced by each upstream signaling feature to downstream responses [87]. Feature-based mutual information is a quantity that measures the amount of information which can be predicted about a response, R, when the measurement of a signaling feature, Xi, is known. The quantity is calculated for each feature to any downstream response at a given time point. Feature based mutual information is a universal quantity by definition, and can be applied to any network topology regardless of its underlying physical basis or complexity. In biological networks, the network topology can be defined as a channel (Fig 1A) or a black box that maps the input functions of signaling dynamics to a downstream response. Because the distribution of the input signals are typically not known a priori, they can be considered as a function of random variables (signaling features) given as X in Eq (9). For a input signal S = fs (X, t), and downstream response R, mutual information I is defined as: I (R; X i) ≡ H (R) - H (R | X i) = H (X i) - H (X i | R) = I (X i; R) (15) where entropy H is H (R) ≡ E p { - log (p (R) ) } = - ∑ R p R (r) log (p R (r) ) and conditional entropy H (R | X i) = E p { - log (p (R | X i) ) } (16) = - ∑ X i p X i (x i) ∑ R p R | X i (r | x i) log (p R | X i (r | x i) ) (17) = - ∑ X i p X i (x i) H (R | X i = x i) (18) Since mutual information is always positive (I (Xi; R) ≥ 0), the conditional entropy is always less than individual entropy, H (R∣Xi) ≤ H (R) [87]. The resulting mutual information I (Xi; R) is then a concave function of pR (r) for fixed pR∣Xi (r∣xi) and a convex function of pR∣Xi (r∣xi) for fixed pR (r). The conditional distribution characterizes these measurements. I (X i; R) = E { log (p X i R (X i, R) p X i (X i) p R (R) ) } (19) = ∑ R ∑ X i P R | X i (r | x i) P X i (x i) log (P R | X i (r | x i) ∑ P R | X i (r | x i) p X i (x i) ) (20) = ∑ X i P X i (x i) (H (p X i (x i) ) - H (P (R | X i) ) ) (21) The uniform distribution of signaling features Xi ensures the maximum possible entropy of H (P X i (x i) ) = log (1 m a x (X i) − m i n (X i) ). The intuition that the conditional distribution P (R∣Xi) is narrower or more concentrated than PXi (xi) is quantified by the fact that the entropy H (P (R∣Xi) ) is smaller than H (PXi (xi) ), and this reduction in entropy is exactly the information when observing R provides about Xi, measured here in bits [88]. However, it should be noted that this equation allows us to measure the information carried by the response to recover estimates of the signaling features in each cells. We then assumed that the total as well as conditional distribution of downstream responses emerging from signaling activation followed a Gaussian distribution, whose variance was not dependent on Xi. In mathematical terms, this means V[R∣Xi] = EXi (V[R∣Xi]). From the law of total variance V[R] = V (E[R∣Xi]) + E (V[R∣Xi]), mutual information is estimated as: I (R; X i) = H (R) - H (R | X i) = 1 2 log (2 π e V [ R ]) - 1 2 log (2 π e V [ R | X i ]) (22) = - 1 2 log E (V [ R | X i ]) V [ R ] = - 1 2 log V [ R ] - V (E [ R | X i ]) V [ R ] (23) = - 1 2 log (1 - V (E [ R | X i ]) V [ R ]) (24) where entropy (H) of any Gaussian distribution pZ (z), is H (Z) ≡ 1 2 log (2 π e V [ Z ]) [60,87]. The information content between R and Xi in Eq (24) was estimated using Eq (13). Since both the signaling feature and and the output response measure are assumed to be a Gaussian, the estimation of the mutual information gives us a lower bound on the information carried by a signaling feature to the downstream response [89]. We later make a generalization that the information is not conveyed by just one single signaling feature, but by a combination of features, as it was the case revealed during the analysis. It is also important to emphasize that the number of bits of information carried by a single feature has a meaning irrespective of the dynamic multiplexing resulting from the network topology. We employed decision tree analysis by generating 105 sets of simultaneous variation in the signaling features through application of Sobol sequence ‘MatousekAffineOwen’ algorithm to scramble a total of 18-dimensional pseudo-random numbers generated using Sobol’s method in MATLAB computational environment [84]. The responses to these variations were statistically grouped into three categories: High, Mid and Low phenotypes. High phenotype represented most positive deviation from the mean of the downstream response at the time point under consideration, Low phenotype meant the most negative deviation from the mean of the downstream response, and Mid represented the median of the response distribution. We then selected first 250 profiles from these groups to identify the rules of the binary tree, using the classification and regression trees (CART) algorithm [43]. The CART algorithm uncovers the predictive structure of how cells might be classified in various phenotypic groups based on variations in dynamical features of the upstream signals. The objective is to progressively split the cells (or predictor dataset) into smaller and smaller subsets until each subset corresponds to only one phenotype (or a leaf of the tree). Each subset represent a hierarchy of binary rules defining the ranges of variation in the signaling features determined through a splitting criteria. In our study, the Gini index was used as the splitting criteria to construct the trees. The analysis was performed in the R statistical language [90] using the ‘rpart’ and the‘rpart. plot’ packages [91,92].
Single cell studies have shown that differential patterns in the dynamics of signaling proteins, transcription factor activity, gene expression, etc. produce distinct downstream outcomes. The opposite also holds true where particular cellular outcomes have been found to be associated with the dynamical pattern of one or more signaling molecules. Signaling pathways, therefore, serve as signal processing units to inform specific downstream regulation. However, the functional capabilities of the dynamic aspects of signaling are not well understood. To address this issue, we developed a new approach that evaluates information processing between dynamic features in signaling patterns and transcriptional regulatory activity. Our work demonstrates that the information transfer occur through decoding of temporal history of signals rather than only through instantaneous correlations. Moreover, our results identify regulatory network motifs as the critical components in the information processing and filtering of variability in signaling dynamics to produce distinct patterns of downstream transcriptional responses. Our methodology can be broadly applied to single cell scale data on experimentally accessible downstream measures to infer dynamic aspects of upstream signaling.
Intracellular Information Processing through Encoding and Decoding of Dynamic Signaling Features