Signal transduction controls heterogeneous NF-κB dynamics and target gene expression through cytokine-specific refractory states

Cells respond dynamically to pulsatile cytokine stimulation. Here we report that single, or well-spaced pulses of TNFα (>100 min apart) give a high probability of NF-κB activation. However, fewer cells respond to shorter pulse intervals (<100 min) suggesting a heterogeneous refractory state. This refractory state is established in the signal transduction network downstream of TNFR and upstream of IKK, and depends on the level of the NF-κB system negative feedback protein A20. If a second pulse within the refractory phase is IL-1β instead of TNFα, all of the cells respond. This suggests a mechanism by which two cytokines can synergistically activate an inflammatory response. Gene expression analyses show strong correlation between the cellular dynamic response and NF-κB-dependent target gene activation. These data suggest that refractory states in the NF-κB system constitute an inherent design motif of the inflammatory response and we suggest that this may avoid harmful homogenous cellular activation.

Shown are (a) representative IB-eGFP and p65-mCherry single cell traces, (b) the cross-correlation function between the corresponding IB-eGFP and p65-mCherry single cell traces, (c) the power spectra calculated for IB-eGFP (in green) and p65-mCherry (in red) traces, (d) the correlation between peak IB-eGFP and p65-mCherry oscillation frequency (in red a fitted linear regression line) (e) the correlation between lag times and p65-mCherry oscillation frequency. 17 C9L cells were analysed (in d and e) and single cell traces were normalized to levels at 0 min, trajectories beyond 120 min were analysed. responses each data pectively, f s to a set from left to right). Cells classified as "responders" base on the IB-eGFP gradient at the time of the 2 nd peak depicted with green colour. Figure 9. IL-1 stimulation induced a refractory state. (a) C9 cells stimulated with IL-1β (depicted with a pink bar) and IL-1β neutralising antibody. Shown is mean and s.d. of 35 single cell traces. (b) C9 cells stimulated with a single 5 min IL-1β pulse (wash performed with media containing IL-1β neutralising antibody at 5 min after stimulation). Shown is mean and s.d. of 37 single cell traces. (c) Responses to two 5 min IL-1 pulses at 60 min interval (washes performed with media containing IL-1β neutralising antibody). Shown is the heat map of individual cell traces (Y axis), clustered by the normalized total IB-eGFP single-cell trajectories (from 39 single cells). Data normalized to the IB-eGFP level at the time of 2 nd pulse, time depicted in minutes. Brokenlines indicate "responding" (above) and "non-responding" (below) cell groups. (d) Responses to a 5 min IL-1 pulse followed by a 5 min TNF pulse at 60 min (first wash performed with media containing IL-1β neutralising antibody). Shown is data from 42 single C9 cells, IL-1 and TNF pulses depicted with pink and blue bars respectively. (e) Fraction of responding cells from data in c and d (shown is mean and data range of two replicates). . FITC-conjugated TNF (top panel, in green) was applied at 0 min for 5 min, washed with media, and fluorescence intensity measured at 15 min after the start of the experiment. Tx-Red-conjugated TNF (middle panel, in red) was applied on the same cells at 60 min, washed with media, and then fluorescence intensity measured at 75 min after the start of the experiment. Corresponding bright field images are shown at the bottom. Scale bar 20 m. Figure 11. Quantification of TNF internalization levels. WT cells were stimulated with two pulses of fluorescently labelled TNF at 60 min interval. Tx-Red-conjugated TNF was applied at 0 min. FITC-conjugated TNF was applied on the same cells at 60 min after the start of the experiment. Shown are the total fluorescence levels in 15 cells measured at 15' and 75' min after the start of the experiment.     Supplementary Figure 24. Global sensitivity analysis of NF-B oscillations in response to continuous TNF stimulation. Shown is the correlation between sensitivity scores for the amplitude (on Y-axis) and the period (on Xaxis) of NF-B nuclear oscillations. Most sensitive parameters (labelled on the graph) include those related to IB and A20 transcription, translation and halflife, NF-B transport, total level of NF-B and IKKK.  Relative sensitivity (RS) performance with respect to nuclear NF-B level for all model parameters (calculated as an average temporal sensitivity coefficient over the 300 min trajectory). Names of most sensitive parameters indicated on the graph.

Supplementary Figure 28. Sensitivity analysis of the TNF and IL-1 receptor modules.
Shown is the comparison of sensitivity coefficient calculated with respect to receptor parameters (total TNF and IL-1 receptors (Rtot), receptor activation rate/cytokine binding rate (r1); receptor de-activation rate (r2), half maximal receptor activation (st)) and IKKK-related parameters (total IKKK in TNF and IL-1 pathway (IKKKtott and IKKKtoti, respectively) and A20 degradation rate (c4). Sensitivity coefficients calculated based on the fraction of 300 simulated cells responding to cytokine pulses at a 60 min interval (as indicated on the graph) with respect to 10% parameter change.
Supplementary Figure 29. Schematics of the TNF equilibrated pulses experiment. Cells were stimulated with 5 min TNF pulses at 70 min interval (p70p) in two sets separated by a 4h equilibration period. Shown is a normalized IB-eGFP trajectory of a representative single C9 cell. Timing of the TNF pulses represented with blue bars.    Responding and non-responding cell simulated with total IKKKt equal to 1.5 x 10 5 and 7 x 10 4 molecules, respectively. S. 28 Deterministic Total IKKKt~N(100000,30000 2 ) Supplementary Table 6. Models and parameters used for specific simulations performed in the study.

Generation of the IB-EGFP BAC
The IB BAC (CTD-3214F11, obtained from Life Technologies, UK) was maxiprepped using Machnery-Nagal Nucleobond BAC100 kit and 100ng was used to transform SW102 E. coli cells (a kind gift from Neal Copeland) by electroporation Colonies were screened by Pulsed Field Gel Electrophoresis for clones that cut with an additional SalI site (inserted between 5' Homology arm and GFP sequence). A single positive clone was detected from ~80 screened. Sequencing confirmed the correct recombination.

Analysis of IB and NF-B p65 levels in BAC cell lines.
The temporal NF-κB response was investigated in developed BAC IκBα-eGFP cell lines. In C9 cells, which were used for most of the single cell experiments, TNFα stimulation resulted in a rapid degradation and re-synthesis of IκBα protein, similar to that in WT cells ( Supplementary Figure 1). At 15 min after treatment, the endogenous IκBα protein was undetectable, which was potentially consequence of low resting levels, while a residual expression of IκBα-eGFP was present. A simultaneous recovery of the endogenous and IκBα-eGFP protein was observed, with a peak at 60 min (after continuous TNFα stimulation). Further accumulation of IκBα levels was observed 120 min after a single TNFα pulse, but not after continuous stimulation, reaching pre-stimulation steady-state (i.e., high IκBα-eGFP and low endogenous IκBα levels). This behaviour was in a good agreement with single-cell data in Fig. 2a.  Figure 2). We expect that this be due to transduction of fluorescently labelled p65-mCHerry, which might enhance production of or stabilise IκBα protein 1,9 .
We found that NF-B p65 oscillation patterns in our stable cell lines were consistent (period of 106 ±30 and 109 ±27 min in C9 and C9L, respectively, as defined by individual peak-to-peak timing of IB-EGFP, Supplementary Often the IB-eGFP fluorescence intensity remains below its initial value as it oscillates (Fig. 1c). We expect that this effect could be due (at least in part) to the time taken for folding and maturation of the fluorophore 10,11 , despite IB-eGFP protein level reaching the pre-stimulated steady state (or higher) at 60 min after treatment (Supplementary Figure 1). Thus, only a proportion of stimulusinduced IB-eGFP molecules might be detected in the microscopy experiment and this would be lower at higher oscillation frequencies. In addition we observe that C9L cell line is characterised by lower amplitude of IB-eGFP oscillations (comparing with C9 cells, Fig. 1c). We expect that this could be due elevated IB-eGFP levels in C9L cells, which could lead to more damped oscillation patterns. Damped oscillations at the single cell level are consistent with partial IB degradation observed in the population (Supplementary Figure 1).

IB and p65 oscillate out of phase.
In response to continuous TNF stimulation C9L cells exhibited out-of-phase oscillations in the IB-eGFP level and the p65-mCherry nuclear translocation

Supplementary Note 2: Single cell responses to TNF and IL-1 pulses
C9 cells were stimulated with two 5 min TNF pulses at intervals ranging from 50 to 100 min (Fig. 2e). The average IB-eGFP traces corresponding to the "responding" and "non-responding" clusters at 60 min interval are shown in Supplementary Figure 6. Single cell IB-eGFP trajectories are depicted with heat maps (see Fig. 2c for 60 min, and Supplementary Figure 7 for 50, 70, 80 and 100 min pulse interval). Cells were then classified as "responders" if a gradient of the corresponding IB-eGFP trajectory at the time of stimulation was not positive; otherwise they were called "non-responders".
Unsupervised k-means clustering was used to validate this classification (Supplementary Figure 8). Cells were first clustered into different number of groups varying from two to four. The previously inferred "responding" and "nonresponding" cells (based on the IB-eGFP gradient) were mapped on to obtained clusters (depicted in green and blue respectively). We found that overall the gradient-based classification was recapitulated very well in the kmeans clusters (with exception of the 50 min pulsing). The Silhouette coefficient could be used to find the optimal number of clusters, however in practice it would be more difficult to implement. For example, 2 optimal clusters were identified in the 60 min pulse interval data (indicated by the highest Silhouette coefficient, in comparison with the coefficients obtained for three and four clusters, respectively), while three of four clusters separated the "responding" and "non-responding" cells better. Based on this analysis, we therefore applied a method based on the IB-eGFP gradient in this manuscript.
In order to test whether multiple refractory states exist within NF-B system C9 cells were stimulated with two pulses of IL-1. 76% of cells responded to the second pulse (Supplementary Figure 9), suggesting existence of a similar stimulus specific refractory state to that seen with TNF. All the cells stimulated with a pulse of TNF responded to the pulse of IL-1 60 min later, suggesting that the refractory period to TNF stimulation might be controlled upstream of IKK in the signal transduction (Fig. 3a-d). When a TNF pulse was applied after IL-1 stimulation, 54% of cells responded, Supplementary Figure 9 (in comparison almost ~100% of cells responded to a IL-1 pulse 60 min after stimulation with TNF, Fig. 3a). This suggested some influence of the IL-1 transduction pathway on TNF signalling.

Supplementary Note 3: Analysis of TNF internalization.
Cells were stimulated with 5 min pulses of FITC and Tx-Red fluorescently labelled TNF at a 60 min interval (Figs 2e-f). Confocal microscopy images of representative WT cell stimulated first with FITC-and then Tx-Red-labelled TNF showed that internalization pattern was not affected (Supplementary  Figure 13). Functionality of the dynasore was confirmed with an acid wash. This suggested that endocytosis did not affect NF-B system responses.

Supplementary Note 4: Correlation between single cell responses and gene expression
Cells were stimulated with three pulses of TNF, all single cells responded to the 1 st and 3rd pulse (at 0 and 100 min, respectively), but were refractory to the 2 nd pulse applied at 50 min. In contrast, when stimulated with IL-1 at 50 min most of the cells exhibited a corresponding nuclear p65-mCherry translocation (Supplementary Figure 32). The comparison between TNF/IL-1/TNF (TIT) and TNF/TNF/TNF (TTT) stimulated genes (as assessed by Nanostring assay, Fig. 6c), showed a panel of 26 genes that was differentially regulated by IL-1 at a time when cells were refractory to TNF (Fig. 6e). In particular, the members of the NF-B signalling system and a number of pro-inflammatory signalling molecules including cytokines and chemokines were differentially regulated (Supplementary Figure 33).
We asked if the differential response to IL-1 and TNF pulsing (Supplementary Figure 33) could be explained by the patterns of nuclear NF-B in single cells. We quantified the level of nuclear p65-mCherry in C9L cells stimulated as in Fig. 6c (for all conditions except of the single IL-1 pulse). In agreement with Fig. 6a

Model development
To theoretically investigate the encoding of pulse TNF stimulation in the NF-B system, the structure 6 and parameters of our previously developed model of the system was used 1,9,14 . The revised model was fitted manually and validated against the following single-cell and population level data (see Supplementary Model parameters were fitted manually based on the parameters that were specifically developed for the SK-N-AS cell line 1 . We noticed that in the fitted model, parameters corresponding to the rate of A20-mediated inhibition of IKK (a reaction assumed in previous models 6,9,17 ) was negligibly small (had no effect on systems behaviour). Therefore we removed this reaction from the model's structure arriving at the final model structure as presented in Fig. 4a.
The developed model consists of 19 ordinary differential equations and 41 parameters describing formation of complexes and their degradation, transport between nucleus and cytoplasm, transcription, translation and regulation of gene activity (see Supplementary Table 2 for description of model variables,   Supplementary Table 3 for equations and Supplementary Table 4 for parameter values). The heterogonous NF-B responses were recapitulated in the model by distributions of parameters related to TNF (IL-1) signal transduction pathway (Fig. 4d).
The effects of stochastic TNF binding and feedback gene transcription were investigated in respective stochastic formulations of the model (see Supplementary Table 5). All simulation conditions used in the study are given in Supplementary Table 6.

Model validation
The model recapitulated single cell NF-B oscillatory responses (period of ~100 min) to continuous TNF and IL-1 stimulation (Supplementary Figure 14). Sensitivity analyses of the TNF and IL-1 receptor modules showed that while related parameters could be relevant (e.g., total number of receptors, Rtot for TNF pulses), the IKKK and A20 levels were more important for single cell refractoriness (Supplementary Figure 28).

Model formulation
The IKK signalling module: The IKK module incorporates the multi-step transduction pathway downstream of TNFR and IL1R. The respective pathways are complex and not fully elucidated 19,20 , however our data suggested that the refractory period was controlled by molecular events above the IKK and below the receptor level (Fig. 3). The simplest biologically relevant structure that could explain this data was used, effectively reducing this complex network of interaction into a single process (Fig. 4a). In particular, the model assumed previously described interaction between the generic IB kinase kinase (IKKK) and the A20 feedback protein 6,9,17,[21][22][23] (Fig. 4B). An alternative mathematical description of the TNF transduction pathway has also been proposed by 3 and included additional molecular events leading in IKKK activation. In agreement with 6,9 , we assumed that IKKK exists in three states: 'neutral' (IKKKn) that can be activated by TNF, 'active' that regulates IKK phosphorylation, and 'inactive' that constitutively converts into a 'neutral' state The total level of IKKK was assumed to be constant per cell, while the A20 inhibited 'active' IKKK, by facilitating the transition to an 'inactive' state and thus limiting the extent of IKK activity. We assumed that TNFR and IL1R transduction pathways utilized different upstream IKKKs (IKKKt and IKKKi, respectively), which independently converged on the IKK 24 . Similarly, we assumed that A20 inhibited both TNFand IL-1-dependent responses. This effect has been described by a prolonged NF-B activation in A20 deficient cells in response to TNF [2][3][4]16 . Signalling downstream from IL-1 is less characterized; in particular no effect was shown by 2,3 in A20 deficient cells, while a prolonged NF-B activation was observed by 16 using knockout studies. For simplicity we used a parallel A20-dependent activation mechanism for IL-1 and TNF transduction pathway.
Distributed IKKK: The distribution of refractory periods was recapitulated in the model by heterogeneous expression of IKKK in the population of cells (Fig. 4d).
Specifically, we assumed that the total IKKK level was constant in a single cell but normally distributed in the population with mean μ=1x10 6 and μ=1.3x10 6 for TNF and IL-1 pathways respectively, and standard deviation =0.3 molecules: , , When simulating the effect of parameters other than IKKK level (e.g. IKKK recycling rate, m3, and IB half life, c4a) on the refractory period, an above normal distribution was used with mean c4a equal to its nominal value and m3 equal to 1.2×nominal value. A standard deviation of =0.3 was used in all cases (as in Fig. 4i). However, other distributions (e.g., log-normal or uniform) could also be used to fit the data (Supplementary Figure 22). to 300 IB and 600 A20 transcript on average in mouse fibroblast cells 6 . In agreement, recent smFISH measurements of A20 and IB transcript abundance in HeLa cells following 10 ng ml -1 TNF stimulation showed heterogeneous expression patterns with average levels of 250 molecules per cell, respectively 28 .
In contrary, our quantitative Nanostring assay suggested that in SK-N-AS cells, the average A20 transcript expression was about 40 times lower than the expression of IB following treatments with TNF and IL-1 (Fig. 6e).
Therefore, while not changing any IB-related parameters, the A20 transcription rate was decreased accordingly, resulting in approximately 6 mRNA and 3,000 protein molecules of A20 on average. Such relatively low levels might generate substantial intrinsic noise (as shown in in the stochastic formulation of the model, Fig. 5b, c). However, the analysis of equilibrated TNF stimulation showed that the NF-B responses were seemingly not affected by this heterogeneity, and instead were deterministic (Fig. 5d). This suggested that the intrinsic noise in A20 regulation has a limited effect on the system's response to TNF pulses. We expect, that the tight control of the enzymatic activity 16,19,20 , rather than apparent heterogeneous expression level drive the function of A20

Stochastic model description
Stochastic simulation algorithm was used to investigate the role of the intrinsic noise in the system. A hybrid model was developed in which the biochemical reactions were split into fast and slow; fast reactions were approximated by rate equations, while slow reactions were considered stochastic 29 . In particular, these former included the regulation of IB and A20 gene activity 9,30 and TNFR receptor binding 6 . In addition, the A20 mRNA death and birth process was considered stochastic due to low number of molecules involved (<10 mRNAs).
We assumed that each gene has two independent homologous copies activated due to NF-B binding and inactivated by the nuclear IB protein 30 . Let denote the binary state of a promoter of single gene copy: whenever NF-κB is bound to the promoter, and G i (t)=0 otherwise. The gene state describes the state of both homologous copies and is given by and for A20 and IκBα promoters, respectively. Let the binding propensity r b be the rate at which the NF-B binds to its unbound promoter, i.e., if the binding site is unbound then the probability of binding in a short time interval of duration dt is r b dt. The dissociation propensity r d is defined similarly. It is assumed that for NF-κB r b is proportional to the amount of nuclear NF-κB, while the dissociation propensity r d is proportional to the amount of nuclear IκBα protein: where G is the state of the gene, and q1x (x={i, a}r for IB and A20 respectively) is the binding rate of NF-B and q2 is the inducible IB dissociation rate.
Because the A20 mRNA transcript was detected at later time points than the IB mRNA in SK-NA-S cells stimulated with TNF 1 , we assumed different NF-B binding rates for A20 and IB promoters, while dissociation rates were kept constant.
When considering an A20 mRNA birth and death process, let rt t and rt d be propensities of synthesis and degradation of one A20 mRNA molecule, respectively. Propensity of synthesis if proportional to the transcription rate and the number of active genes, while propensity of degradation is a function of number of A20 mRNA molecules (tA20): • , Finally, the TNFR receptor-biding propensity, rr b is proportional to the number of available TNF molecules. Similarly, receptor inactivation propensity (due to dissociation of TNF), rr d is proportional to the number of active receptors (Rt): • , • .
The total propensity function r(t) for the occurrence of any of the reactions is given by where rA20 and rIκBα are individual propensities for A20 and IκBα, respectively.
The numerical scheme implemented follows that of 30 : 1. At simulation time t, for given state of the system calculate the total propensity function r(t) of occurrence any of the reactions. where ri(t + τ), i=1,..., are individual reaction propensities and k is the index of the reaction. Time t+τ is then replaced by τ and the process is repeated.

Model simulations
Deterministic model: The system was initialized with the total NF-κB in the complex with IκB in the cytoplasm, and the total IKK and IKKK was initialized in the respective neutral form (all other state variables are set to zero). The total IKKK level was randomized from the reference distribution per cell, and the log( p1)  r(s)ds  0.
system was simulated to reach the steady state before TNF (and IL-1) activation.
Stochastic model: The stochastic receptor binding was considered to recapitulate low dose TNF simulation. Separately, the intrinsic noise due to activation of A20 and IB genes we preformed to show that such model could recapitulate TNF pulse data at 70 min interval (Fig. 5b). This model, however subsequently failed to agree with the equilibrated TNF pulse data (Fig. 5c-e), suggesting that extrinsic noise was more important. In the intrinsic noise simulations, cells were simulated for time t uniformly distributed on the interval from 2 to 5 hours prior to TNFα stimulation. Due to the natural degradation of IκB and the resulting basal NF-κB translocations sensitivity to any initial conditions is lost 30 .
In the stochastic receptor activation model, reactions in equation 14 (Supplementary Table 3) were replaced by propensities for receptor binding and dissociation, and simulated using the numerical scheme described above. In this description, active receptor number (Rt) became a discrete random variable.
In the stochastic transcription model, reactions in Equation 10 (Supplementary Table 3) were replaced by propensities for A20 gene activation and mRNA transcription. Additionally, Equation 7 refers to discrete state of IB gene (GIB), while equation 11 refers to discrete number of A20 mRNA transcript (tA20): 7. ddt(tIB) = c1a*GIB -c3a*tIB 11. ddt(A20)= c2*tA20 -c4*A20 The resulting system was simulated using the numerical scheme described above (with all other equation as in Supplementary Table 3).

Dynamic sensitivity analysis
Dynamic sensitivity analysis was performed using a method outlined in 18 .
Absolute sensitivity equations were calculated by taking the partial derivative of variables Xi, i=1,..,n with respect to the parameter of interest, θj where k is the time instance and N is the number of points in the output.