Eukaryotic cells sense a wide range of signals via surface receptors and in turn mount specific responses by activating gene expression programs. Although the majority of biochemical information on cell signaling has been obtained from population level studies, it is not clear if population data faithfully reflect how individual cells respond.2, 3
For example, pulsed responses of p53 to radiation damage are evident only at the single cell level, and are blurred out in population measurements.4
Similarly, a recent single cell study of LPS induced NF-κB signaling with showed that only half of cells responded to the secondary TNF-α. autocrine signal, creating distinct subpopulations.5, 6
Determining the variation at the single-cell level is turning into a powerful tool both for understanding drug response7
and for general understanding of how biological systems work. Data obtained from single-cell culture measurements is often complementary to what one obtains from in-vivo imaging where cells exist in more complex contexts such as three-dimensional contact and tissue specific signaling environment.
To investigate how individual cells respond to variation in input signal level, we studied nuclear localization dynamics of the transcription factor NF-κB under stimulation with the inflammatory signaling molecule TNF-α. NF-κB is a highly conserved component of the eukaryotic immune system.8, 9
It controls the expression of hundreds of genes in response to a wide range of stimuli including physical stress, UV light exposure, signaling molecules, and pathogens such as bacteria and virus. The dysregulation of NF-κB is involved in a variety of pathologies such as chronic infection, cancer, inflammatory disease, autoimmune disease and improper immune system development. Population studies10
have not revealed the intricate network of information one observes at the single cell level.11, 12
Previous single cell experiments were limited to high TNF-α concentrations (10 ng/ml) and relatively small statistics (~200 cells). NF-κB pathway activation and dynamic properties at lower doses have remained an open question.
We studied 3T3 mouse fibroblast cells6
in a microfluidic cell culture platform13
to measure NF-κB activity under 10 different TNF-α concentrations (100 ng/ml to 0.005 ng/ml) with single-cell resolution (Supplementary Fig. S1
). When the cells are stimulated, NF-κB is transported from the cytoplasm to the nucleus and back out again in characteristic oscillations which we observed with a fluorescent fusion protein (See Supplementary Movies 1, 2
). The microfluidic system allows measurement of all 10 TNF-α concentrations side-by-side in a single experiment with excellent reproducibility, and mimics physiological conditions in terms of volume, fluid flow and number of ligands more plausibly than conventional culture environments where secreted signaling molecules are quickly diluted into milliliters of media. More than 400 live cells were quantified at each condition () and each experiment was repeated four times, extending the throughput of previous such measurements by more than an order of magnitude (Supplementary Table 1
). In parallel experiments with conventional (population based) culture, we used rt-PCR and digital-PCR to quantify time dependent expression profiles of 23 target genes under the same concentration range to link NF-κB dynamics to transcription.
Fig 1 A Real-time fluorescent images of live cells during stimulation with 10 ng/ml (top row) and 0.25 ng/ml (bottom row) TNF-α. Arrows show the activated cells. At the high dose, all cells except one respond, while only two out of five respond at the (more ...)
Our measurements reveal the response characteristics of cells at different signal intensities (i.e. external TNF-α concentration). One of the most intriguing findings is the discrete nature of single cell activation: not all cells responded to TNF-α (, Supplementary Figure S2
and Supplementary Movies
), and the fraction of activated cells decreased with decreasing TNF-α dose (). While nearly all of the cells are activated at doses above 0.5 ng/ml, the percentage of responding cells fell to 50% for 0.1ng/ml and fell to 3% at 0.01ng/ml. The mean fraction of activated cells of four experiments (N~5000 single cells) is well described by a Hill function with n=1.5. The activated single cell traces can be seen in . Strikingly, the median nuclear intensity changed by only a factor of four in response to a four order of magnitude reduction of the TNF-α concentration ( and Supplementary Figure S3
). Similarly, the mean area under the first peak remained relatively constant over the entire range of concentrations, showing that total NF-KB nuclear activity was equal during the primary response (), indicating digital activation in response to TNF-α.
Fig 2 NF-κB response characteristics at different TNF-α doses. Error bars indicate the standard deviation from the median response time, and are indicative of the natural variation in single cell response. A Peak nuclear p65/dsRed localization (more ...)
To determine whether the heterogeneous cell activation is entirely based on a pre-existing variation in cell sensitivities14, 15
(so-called extrinsic noise16
), we stimulated our cells with two 20 minute pulses of intermediate dose (0.1 ng/ml) TNF-α, while allowing the cells to reset by giving a 170 minute break between pulses. We found that 10 % of the cells respond to only the first pulse, 9 % respond to only the second pulse, and 11 % respond to both pulses (Supplementary Figure S4
). The existence of cells responding to only one of the pulses shows that the activation is governed partly by a stochastic element. However, the probability of responding to both pulses is highly increased (~10 fold) compared to a purely stochastic process, which shows that pre-existing variation in internal variables of the cells also plays an important role in the outcome.14, 15
The heterogeneous activation of single cells requires the re-interpretation of previous population level studies and the models derived from them.10
The population average of only active cells yields an accurate mean nuclear NF-κB intensity (, bottom right). On the other hand, the average over the entire population including active and inactive cells shows significantly reduced mean intensity and completely washed out dynamics, highlighting the importance of single cell measurements.
Although the cells have a stereotyped digital response when activated, there are aspects of their response that are dose dependent and therefore represent analog information processing of the stimulus. Most importantly, the cells show significantly delayed activation under lower doses (). The activation response time (the time from TNF-α stimulation to the first peak) is around 50 minutes for lowest doses, while the cells activate within 20 minutes at highest doses around 10 ng/ml. The variation in response time at the highest doses is close to the experimental limit of 6 minutes, and increases to about 60 minutes at the lower doses. In contrast to the response time, the mean decay time of the first nuclear localization peak is surprisingly independent of the TNF-α dose, around 20 minutes for all concentrations (Supplementary Figure S5
A related dose-dependent parameter is the number of nuclear oscillations, which also depends on the initial TNF-α concentration (). The intensity of the secondary NF-κB peaks remains relatively constant at high and intermediate doses (Supplementary Fig S6
), but they completely disappear in some cells at the lowest doses measured. The time between the first and second peaks remained constant at 80 minutes over the measured range (Supplementary Fig S7
). Later peaks quickly lose synchrony, and the significant de-phasing after the third peak results in washing-out of the single cell dynamics at the population level (Supplementary Fig S3
). These subtle variations in the dose-dependent response provide strong constraints on mathematical models of TNF-α induced NF- κB signaling, which are discussed in more detail below.
At saturating doses of TNF-α, it is known that NF-κB nuclear oscillations drive expression of three classes of genes: early, intermediate and late.11, 12, 17–19
We measured the time dependent expression levels of 23 NF-κB target genes on populations of cells from these classes under a broad range of TNF-α doses using quantitative rt-PCR and digital-PCR ( and Supplementary Figure S8
). A series of genes are up-regulated at early times, including IκB-α and A20 (). Transcripts of early genes have high degradation rates,19
which results in expression profiles that closely follow the NF-κB nuclear localization dynamics (). An intriguing property of the early genes is that the expression levels normalized to the fraction of active cells in each population is relatively constant for all doses (), while it is strongly dependent on dose for late genes (). This suggests that the initial NF-κB localization burst is sufficient to express fast genes even at the lowest stimulation doses, indicating digital signaling to early gene expression in agreement with digital activation of NF-κB shown in .
Figure 3 Time dependent expression profiles of NF-κB target genes under TNF-α stimulation with doses ranging from 10ng/ml to 0.01ng/ml, measured on cell populations containing nearly 500 cells in each measurement. Relative expression levels were (more ...)
The expression levels of the intermediate and late genes build up slowly following persistent NF-κB oscillations (). Such oscillations require constant TNF-α stimulation and thus are observed only when initial TNF-α level is highest (). The initial localization peak seen in all active cells is not sufficient to induce the late genes. Among the late genes that are highly expressed, CCL2, TRAIL, CASP4, RIPK2 and RANTES lead to apoptotic signaling and Serpina 3g acts to prevent apoptosis.20
This suggests that cells that are continuously exposed to high dose TNF-α activate pro and anti-apoptotic gene expression programs mediated by late oscillations of NF-κB, while low dose TNF-α signals do not lead to apoptosis related signaling due to lack of late oscillations. Furthermore, other signaling molecules and regulatory proteins are also up-regulated with TNF-α stimulation, including IL-6, CCL-20, ICAM-1, and ZFP36, suggesting possible synergistic or interfering mechanisms.
The wealth of biochemical data has resulted in a number of mathematical models that describe TNF-α induced NF-κB oscillations. Hoffmann et al. created a model of the core NF-κB network with negative feedback from three known IκB isoforms.18
Nelson et al. modified and fitted the same model to single cell data to show that NF-κB oscillations play role in regulating gene expression.11
However, this model is deterministic and is unable to reproduce the heterogeneity in single cells; it also predicts extremely high IκB transcription rates that are not biologically plausible. Lipniacki et al. modified the model by incorporating the effect of A20 which negatively regulates NF-κB; 21, 22
and using stochastic dynamics for slow reactions.23, 24
This model is simplified in that it omits negative feedback from IκB-β and IκB-ε, yet it succeeds in reproducing NF-κB activity at the single cell level with the natural heterogeneity, and uses biologically plausible internal parameters. A version of this model with modified A20-IKK interaction structure was able to explain the response of the network to repeating short pulses of TNF-α in SK-N-AS cells.12
The bimodal nature of single cell response we observed (i.e. whether the cell activates or not) and the overall response characteristics shown in are not simultaneously present in the existing models. We therefore built a mathematical model () to explain our observations based on the model architecture by Lipniacki et al.24
The original model assumes stochastic ligand-receptor and transcription factor-gene binding, which was preserved in our model, and new constraints that stemmed from our data were introduced. Calculations of ligand-receptor binding at low concentrations in the microfluidic environment show that TNF-α trimers are quickly bound by thousands of receptors on each cell (Supplementary Fig. S9
). This results in the rapid loss of extra-cellular TNF-α and a lower number of active cells in the microfluidic chambers compared to culture wells. This effect was introduced to the model as concentration dependent loss of extra-cellular TNF-α through receptor internalization. Furthermore, we have included variation in TNF-α receptor levels to account for the variable sensitivities we observed in short-pulsed stimulation of our cells as shown in Supplementary Figure S4
Fig 4 A New mathematical model architecture based on stochastic receptor and gene binding and second order nonlinear IKK activity. B Model simulations of NF-κB nuclear localization at the experimentally measured TNF-α concentration range. The (more ...)
Nevertheless, the combination of stochastic receptor binding and ligand loss were not sufficient to produce switch-like activation observed at low doses. The linear model architecture results in the gradual increase in the nuclear intensity with increasing TNF-α dose, and failed to produce a clear separation between active and non-active cells. We assumed a nonlinear IKK (inhibitor of NF-κB kinase) activation profile, where the activation rate of the IKK depends quadratically on the concentration of the upstream kinase IKKK. Such a nonlinear activation rate is justified by the fact that IKK subunits IKK-α and IKK-β must be phosphorylated at two serine residues, S177 and S181 to attain full activity.25
Incorporation of the quadratic activation resulted in achieving the bimodal switching similar to our experiments.
The model () was manually fitted to our data, and a single set of biochemically measured, assumed and fitted parameters were used for all the simulations (See Supplementary Mathematical Methods
). The core model comprises 16 ordinary differential equations representing biochemical reactions and 34 rate constants, where 20 of the constants are fixed (biochemically measured or based on related data), and the other 14 were varied for model fitting within biologically reasonable boundaries. We choose the set of parameters which not only fits to our data, but also produces qualitative agreement with a major subset of the existing literature data. The internal parameters of the model, such as the number of proteins and mRNAs, production, degradation and transport rates are all biologically plausible values.
The simulated single cell traces for 3T3 cells at different TNF-α doses can be seen in , which agree well with the data in terms of timing, intensity and the variability of the first peaks at all doses. The number and timing of the secondary peaks are reproduced at the highest dose with the current fitting parameters; however at lower doses the current model parameters do not successfully reproduce the number and timing of secondary peaks. Most remarkably, the model can reproduce -within a single set of biochemical parameters- the measured concentration dependent single cell activation probability (), mean nuclear NF-κB intensities () and the response times () and their distributions. Furthermore, it reproduces with excellent accuracy the fraction of cells responding to consecutive short pulses of low dose TNF-α (Supplementary Figure S12
and Math. Methods Table 3
). The extrinsic variability in NF-κB protein and TNF-receptor levels increases the heterogeneity in nuclear intensities and peak timings by approximately two-fold compared to simulations with only intrinsic noise, but does not change the mean. The combination of extrinsic and intrinsic noise result in a better fit to the experimental variation we observed.
Simulated expression profiles for early, intermediate and late genes can be seen in . The gene expression closely follows NF-κB dynamics and agrees well with most of the genes we measured in . Some measured expression profiles such as IL 6, TRAIL and CCL2 do not follow NF-κB localization strictly and the single transcription factor model we used was not able to reproduce these profiles, which suggests that parallel synergistic or interfering pathways are in effect in transcription of these genes.27
In conclusion, we have mapped the localization dynamics of NF-κB under persistent TNF-α stimulation by measuring thousands of single cell traces, and time dependent gene expression profiles over a concentration range of 4 orders of magnitude. A mathematical model that reproduces NF-κB dynamics and transcription at the measured concentration range was developed, constituting a significantly improved model of the NF-κB pathway under TNF-α. We have found that cells activate in response to TNF-α in a switch-like manner. Although the fraction of activated cells decreases to zero with decreasing TNF-α dose, the NF-κB amplitude remains high, allowing high expression of early genes. Early gene expression is not dependent on the inducing signal intensity, while the expression of later genes requires persistent nuclear localizations of the transcription factor, which exist only at the highest input signal levels. The cells process analog signal intensity information using a diverse set of parameters such as the timing, intensity and number of transcription factor oscillations to create digital outputs (activation and early gene expression). In addition to their biological significance, our findings highlight the importance of high-quality, single-cell data in understanding and modeling biological systems, and demonstrate the efficiency of high-throughput microfluidic techniques in obtaining such data.