|Home | About | Journals | Submit | Contact Us | Français|
Cells operate in ever changing environments using extraordinary communication capabilities that emerge from the interactions of genetic circuitry. The mammalian immune response is a striking example of the coordination of different cell types.1 Cell-to-cell communication is primarily mediated by signaling molecules that form spatiotemporal concentration gradients, requiring cells to respond to a wide range of signal intensities. Here we use high-throughput microfluidic cell culture, quantitative gene expression analysis and mathematical modeling to investigate how single mammalian cells respond to different concentrations of the signaling molecule TNF-α, and relay information to the gene expression programs via the transcription factor NF-κB. We measured NF-κB activity in thousands of live cells under TNF-α doses covering four orders of magnitude. We find, in contrast to population studies, that the activation is heterogeneous and is a digital process at the single cell level with fewer cells responding at lower doses. Cells also encode a subtle set of analog parameters to modulate the outcome; these parameters include NF-κB peak intensity, response time and number of oscillations. We developed a mathematical model that reproduces both the digital and analog dynamics as well as the most gene expression profiles at all measured conditions, constituting a broadly applicable model for TNF-α induced NF-κB signaling in various types of cells. These results highlight the value of high-throughput quantitative measurements at the single-cell level in understanding how biological systems operate.
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 (Fig 1) 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.
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-α (Fig 1A, Supplementary Figure S2 and Supplementary Movies), and the fraction of activated cells decreased with decreasing TNF-α dose (Fig. 1B). 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 Figure 1C. 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 (Fig 2A 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 (Fig 2B), indicating digital activation in response to TNF-α.
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 (Fig 1C, 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 (Fig 2C). 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 (Fig 2D). 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 (Figure 3 and Supplementary Figure S8). A series of genes are up-regulated at early times, including IκB-α and A20 (Fig 3A). Transcripts of early genes have high degradation rates,19 which results in expression profiles that closely follow the NF-κB nuclear localization dynamics (Fig 3B and 3C). 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 (Fig 3D), while it is strongly dependent on dose for late genes (Fig 3E). 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 2.
The expression levels of the intermediate and late genes build up slowly following persistent NF-κB oscillations (Fig 3B). Such oscillations require constant TNF-α stimulation and thus are observed only when initial TNF-α level is highest (Fig 3E). 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 Figure 2 are not simultaneously present in the existing models. We therefore built a mathematical model (Fig 4A) 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.
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 (Fig 4A) 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 Figures 4B, 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 (Fig 4C), mean nuclear NF-κB intensities (Fig 4D) and the response times (Fig 4E) 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 Fig 4F. The gene expression closely follows NF-κB dynamics and agrees well with most of the genes we measured in Fig 3. 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.
We previously used a lentiviral system to create p65−/− mouse fibroblast (3T3) cells expressing the fluorescent fusion protein p65-DsRed under control of the endogenous mouse p65 promoter. 6 Overexpression of p65 can cause unusual NF-KB activation, we therefore cloned the 1.5 kb upstream of the relA gene into the lentiviral construct to control p65-FP expression, and used the resulting construct to infect relA knockout 3T3 cells. These cells exhibited a normal response to TNF at the population level. To relieve a bottleneck in image processing in this study, we also infected the cells with a lentivirus containing the nuclear marker H2B-GFP driven by the human Ubiquitin C promoter. After cloning, the cells were frozen and newly thawed cells were used for each experiment to prevent 3T3 cells transforming and to minimize heterogeneity. A correlation between p65-DsRed levels and cell activation was not observed.
Cells were seeded at constant density at 20,000 cells/cm2 (200 ± 25 cells/chamber) into microfluidic chambers and were cultured for one day to reach near 50 % confluency prior to stimulation experiments. The external conditions were set to standard culture conditions, with 5% CO2 and 37°C external temperature and maintained at this level. During cell growth, 33% of the chamber volume was replaced with fresh media (DMEM) every hour using the nanoliter microfluidic pump, which resulted in vigorous proliferation of 3T3 cells. We made sure that cells were healthy, motile and proliferating before TNF-α stimulation experiments. We were able to establish such cultures for durations up to a week in the microfluidic chambers, depending on initial seeding density. During stimulation experiments, various TNF-α (Roche) dilutions in DMEM were kept at ice at all times and the fluids were immediately flowed into the microfluidic chambers during stimulation to avoid degradation of TNF-α. The high surface-area to volume ratio of the plastic tubing and the microfluidic environment we used allow rapid warming of the media when it enters the environmental chamber, and therefore the media flown in the chambers were at ambient (37°C) temperature. The media and TNF-α vials were pressurized with 5% CO2 to minimize Ph change. The cell viability was not affected negatively during media flow and stimulation. The microfluidic culture chip allows the replacement of the entire culture media in less than a second, which results in a perfect step-like TNF-α concentration profile on the time scales of NF-κB dynamics. Once the TNF-α loaded media was flowed in, the chambers were sealed and were imaged every 6 minutes in both GFP and dsRed fluorescence channels during the entire experiment, for up to 12 hours. The cells were not fed with media after stimulation and during imaging; they remained in the same TNF-α loaded media during the entire experiment without disturbance. All TNF-α concentrations (up to 10 at a time) were tested in parallel chambers (with duplicates) in a single experiment, and such experiments were repeated multiple times with minimal variation between results. We therefore tested up to 20 different microfluidic chambers at a time. To avoid adsorption of TNF into native PDMS, we have used each experimental chamber only once during stimulation experiments. Furthermore, the microfluidic chambers were pre-coated with fibronectin which minimizes adsorption of other molecules.
Cells were imaged in both GFP and DsRed channels with a Leica DMI6000B microscope (20x air objective) and Retiga-SRV CCD camera (Q Imaging) every 5–6 minutes for several hours. Cell traces were visually checked and compared to images to exclude mitotic cells, as mitotic cells typically show increased signal intensity that could be confused with active cells. The fraction of responding cells at each concentration was calculated as the number of responding cells divided by the sum of responding and non-responding cells. The fraction of cells activated and the individual classifications for cells showed a high correlation between two independent measurements. To generate time courses of NF-κB localization in single cells, custom Matlab software using the Image Processing Toolbox was used to automatically identify the nuclei from H2B-GFP images and extract nuclear intensities from p65-DsRed images. To identify nuclear regions, H2B-GFP images were local range contrast filtered with a neighborhood of 3 pixels and then thresholded, where the threshold level was automatically determined by k-means clustering pixel intensities with k = 3. Touching or merged nuclei (determined by solidity < .925) were then separated by a watershed transform with markers seeded at k-means clustered centroids. Nuclei between time points were then linked to the nearest nuclei in the next time point and preliminary quality control checked for constant nuclear area through tracking. H2B-GFP and p65-DsRed images were aligned by cross-correlation when necessary and nuclear intensities were then extracted from the associated p65-DsRed image. For the first time point, cytoplasmic areas and p65-DsRed intensities were extracted through a combination of local thresholding and watershed transforms with the nuclei as marker seeds. All cell tracking was manually checked to eliminate mistakes by the automated analysis. Cells that divided or left the field of view during the experiment were not included in the analysis of NF-κB localization dynamics. Time courses of NF-κB localization were linked to the responding or non-responding classification using the XY coordinates in the field of view. The mspeaks function in the Bioinformatics Toolbox of Matlab was used to identify peaks in NF-κB nuclear localization and to calculate the times of full width half max (FWHM) of the first peak. Localization dynamics from cells that were judged to be non-responders were used to determine a peak threshold for responding cells. Because we could not always image the ascending portion of the first peak with sufficient temporal resolution, we calculated the difference between the time of the first peak (response time) and the later FWHM time, which we called the first peak decay time. Subsequent statistical analysis was performed with Matlab.
The peak intensity and timings of the peaks were found using custom peak detection software written in Matlab. The area under each nuclear localization peak is a measure of the total NF-κB nuclear activity that lead to the transcription of target genes (i.e. the total amount and duration of NF-κB presence in the nucleus), while the peak intensity is a measure of the maximum level of the transcriptional activity. At each dose, the area of the first peak was integrated from the time of TNF-α stimulation to the first minimum for each cell using Matlab. The delayed rise time observed during low-dose stimulation results in increased peak width, which compensates the reduced peak intensity, resulting in nearly equal median first peak areas at all doses.
Cells cultured in well plates were stimulated with 6 different TNF-α concentrations from 0.01ng/ml to 10ng/ml prepared in DMEM media. At the end of each stimulation experiment (0.25, 0.5, 2, 4, 6, 8, 10 and 12 hours after stimulation), cells were lysed and cDNA was synthesized using Cells Direct One Step RT-PCR kit (Invitrogen), and TaqMan primers and probes (Applied Biosystems) were used for real-time RT-PCR. These were population based measurements (we had nearly 500 cells in each RT-PCR reaction per condition), and they contained both active and non-active cells. Each experimental condition was duplicated, and was compared to control wells where cells were not stimulated with TNF-α. Each condition was measured at least three times. We have used microfluidic dynamic arrays (Fluidigm 48×48 Dynamic Array) to perform as many RT-PCR reactions in parallel as possible, and all 64×24 experiments (x4 replicates of each reaction) were performed with only three dynamic array PCR experiments. The cycle thresholds (CT) (Supplementary Table 2) measured during RT-PCR reactions were converted into relative expression levels (2−CT). The relative gene expression levels were then calibrated for total mRNA molecules per cell via digital-PCR measurements performed on a single gene (Fluidigm 12 Digital PCR Chip). More information about the use of digital-PCR for can be found in Anal. Bioanal. Chem. 394, 457–467. Figure 3A, 3B and 3C shows mRNA levels measured on populations that contain active and non-active cells. In Figures 3D and E, we estimate the mRNA levels for only active cells, by dividing the expression levels with the fraction of active cells at each dose measured in Figure 1B. The Figures 3D shows that early genes are expressed at very high levels even at the lowest TNF-α doses, demonstrating how single cell data can be used to interpret population based measurements of gene expression.
See Supplementary Material for a full discussion of mathematical modeling and fitting.
We thank A. Leyrat and R. Gomez-Sjoberg for development of the automated cell culture system, assistance in modifying the software for these experiments, and for contributions to experimental design and preliminary data. This research was supported in part by an NIH Director’s Pioneer Award (to SRQ), an NCI Pathway to Independence Award (K99CA125994, to MWC) and Foundation for Polish Science (TEAM 2009-3/6) (to TL).