Single TNFα trimers mediating NF-κB activation: stochastic robustness of NF-κB signaling
© Lipniacki et al. 2007
Received: 23 February 2007
Accepted: 09 October 2007
Published: 09 October 2007
Skip to main content
© Lipniacki et al. 2007
Received: 23 February 2007
Accepted: 09 October 2007
Published: 09 October 2007
The NF-κB regulatory network controls innate immune response by transducing variety of pathogen-derived and cytokine stimuli into well defined single-cell gene regulatory events.
We analyze the network by means of the model combining a deterministic description for molecular species with large cellular concentrations with two classes of stochastic switches: cell-surface receptor activation by TNFα ligand, and IκBα and A20 genes activation by NF-κB molecules. Both stochastic switches are associated with amplification pathways capable of translating single molecular events into tens of thousands of synthesized or degraded proteins. Here, we show that at a low TNFα dose only a fraction of cells are activated, but in these activated cells the amplification mechanisms assure that the amplitude of NF-κB nuclear translocation remains above a threshold. Similarly, the lower nuclear NF-κB concentration only reduces the probability of gene activation, but does not reduce gene expression of those responding.
These two effects provide a particular stochastic robustness in cell regulation, allowing cells to respond differently to the same stimuli, but causing their individual responses to be unequivocal. Both effects are likely to be crucial in the early immune response: Diversity in cell responses causes that the tissue defense is harder to overcome by relatively simple programs coded in viruses and other pathogens. The more focused single-cell responses help cells to choose their individual fates such as apoptosis or proliferation. The model supports the hypothesis that binding of single TNFα ligands is sufficient to induce massive NF-κB translocation and activation of NF-κB dependent genes.
Living cells are considered noisy or stochastic biochemical reactors. Most of the cell to cell variability is due to existence of stochastic switches or slow reaction channels involving limited numbers of reacting molecules. Stochastic switches provide inputs for amplification cascades, which translate the single molecule events into a larger population of downstream effector molecules.
The most studied example of stochastic regulation is gene expression, where stochasticity, in eukaryotic organisms, arises mostly from fluctuation in gene activity [1–5] and mRNA synthesis or decay [6–10] reviewed recently in . Control of gene activity is mediated by transcription factors that bind to specific promoter regions, switching the gene on or off. When the gene is active, RNA polymerase may bind to the gene promoter and enter the transcriptional elongation mode, producing full length pre-mRNA transcripts. The edited mRNA is then exported to the cytoplasm, where the protein translation occurs. In this way, a single gene activation event results (if the activation period is sufficiently long) in a burst of mRNA molecules , which is then translated into an even larger burst of proteins.
Another example of stochastic regulation is provided by cell surface receptors and amplification of successive cascade of downstream kinases. There is a large body of evidence that cells are capable of responding to a single- or a very limited number of activating molecules. For example, retinal Rod cells are able to transduce a single photon into a hyperpolarization response . For the present study, signal amplification from a small number of receptors detecting pathogen presence is specially important since it enables cells to respond with protective cytokine cascades to protect the tissue before adaptive immune response can be generated. For example, T-lymphocytes are able to detect a single foreign peptide antigen and only three peptides are required for induction of T-cell cytotoxicity [14, 15]. Similar behavior is seen for the toll-interleukin (TIR) superfamily of receptors, where IL-1 signal transduction has been observed in cells expressing only about 10 IL-1 receptors per cell . Similarly, cells respond to TNFα stimulation at femtomolar concentration, i.e., when the number of TNFα ligands per cell is very limited . In all cases, the activation of receptors leads to the initiation of a signal transduction-amplification cascade, involving activation (phosphorylation) of downstream effector kinases. More generally, it is important to point out that cells can detect and respond to single molecule intracellular events such as DNA damage (leading to p53 activation) or presence of single viral RNAs. Although the phenomenon of cell sensitivity to single activating molecules is well established in biological systems, very few theoretical studies have addressed the effect of stochastic cell surface signaling and its consequences for the downstream cellular responses.
Innate immunity is an intensively studied cellular signaling response to pathogens and pathogen-associated patterns which results in the expression of protective cytokines, such as interferon β, that serve to limit the spread of infection until more specific adaptive immunity can be generated. In this regard, the cytoplasmic transcription factor NF-κB is a major mediator of innate immune responses [18, 19]. In resting cells NF-κB is sequestered in the cytoplasm by dimerization with inhibitory proteins called IκB. Although several IκB isoforms have been identified, the primary inhibitor is IκBα. In the classical NF-κB activation pathway, extracellular signals such as tumor necrosis factor-alpha (TNFα) and interleukin-1 (IL-1) bind to cell surface receptors coupled to the cytoplasmic IκB kinase (IKK), a multiprotein complex that phosphorylates IκBα, leading to its ubiquitination and then to its rapid proteasomal degradation . Liberated NF-κB is then rapidly translocated into the nucleus to bind to high affinity sites in the genome, thereby influencing target gene expression. Experimental findings have shown that NF-κB nuclear residence is transient and dynamic, an observation that has led to the discovery of negative feedback NF-κB-IκBα loop in the NF-κB pathway .
The two levels of autoregulatory negative feedback control, termed the NF-κB-IκBα and NF-κB-A20-IKK feedback loops, arise because both IκBα and A20 genes are directly regulated by NF-κB binding. In the NF-κB-IκBα feedback loop, NF-κB enters the nucleus after IκBα degradation, NF-κB induces IκBα resynthesis to recapture activated nuclear NF-κB and return it to the cytoplasm. In the NF-κB-A20-IKK feedback loop, A20, a protein that is not expressed prior to stimulation is also strongly NF-κB responsive . A20 is an inhibitor of IKK kinase that complexes with an IKK regulatory subunit leading to its inactivation and that induces degradation of RIP a necessary component of the active TNFR1 receptor complex. Without A20 expression, the IKK retains activity which leads to rapid degradation of the newly resynthesized IκBα, which destroys the NF-κB-IκBα feedback loop . As an illustration, genetically A20 deficient mice are hypersensitive to TNFα and develop severe in inflammation even though they have an intact IκBα mRNA expression .
Nuclear NF-κB activates groups of genes through a process initiated by its binding to high affinity DNA binding sites in regulatory regions of their promoters. Although NF-κB binding to some genes results in rate-limiting complex formation of coactivators, pre-initiation factors, and RNA polymerase (Pol) II, the mode of regulation of rapidly induced negative feedback inhibitors appears to be distinct. Interestingly, chromatin immunoprecipitation assays have shown that the IκBα and A20 promoters are "pre-loaded"- already bound by general transcription factors, coactivators and RNA Pol II, waiting for the presence of NF-κB binding to activate expression [25, 26]. In these promoters, RNA Pol II is stalled in an activated state; upon NF-κB binding, RNA Pol II enters a competent transcriptional elongation mode and becomes able to transcribe the gene. In this manner, the inhibitory genes of the NF-κB feedback loop are poised to rapidly respond to the presence of nuclear NF-κB.
In last five years several models of the NF-κB signaling regulatory module have been developed, reviewed in . The first attempt was a one-feedback loop model that concentrated on the interplay between the three IκB isoforms . The next attempt was our two feedback loop model , incorporating effects of both the IκBα and A20 inhibitors. In this model the representation of the NF-κB-IκB regulatory module was simplified by incorporation of only one IκB inhibitor, IκBα, that is responsible for the majority of cytoplasmic NF-κB binding and the only one which knockout is lethal, . Incorporating the second NF-κB-A20-IKK negative feedback loop, accurately predicted the time dependent profile of IKK activity. A third model introduced the transduction pathway, from the TNFR1 receptor to activation of the IKKK and IKK kinases and was used to analyze responses of HepG2 cells and HepG2.2.15 (HepG2 cells producing hepatitis B virus) to TNFα stimulation .
Recently, the NF-κB signaling in single cells was analyzed both experimentally [31–34] and numerically by means of stochastic modeling [35, 36], or agent-based modeling . Both analyses indicate that the NF-κB regulatory module demands single cell, stochastic analysis due to cellular heterogeneity and population asynchrony. In the present work we expand our two-feedback loop, single-cell stochastic model of NF-κB activation to incorporate a second stochastic switch at the level of the TNFα-TNFR1 interaction. We analyze response of the NF-κB regulatory module over a broad range of stimulation by its activating ligand. We will show, how, although this may seem counter-intuitive, stochasticity and stochastic switches may introduce robustness into the gene regulatory response. In short, the stochastic robustness is due to amplification cascades and progressive signal saturation. As a result, a low amplitude signal (small concentration of activating ligand or transcription factor) leads to an "almost" Yes or No response, with the probability of Yes being a function of the input signal amplitude. This type of regulation enables cells to chose a well-defined fate such as signaling, apoptosis, or others, but in addition it allows individual cell responses to vary. In the range of stimulation amplitudes, for which most cells follow the same evolution path, the cell population-based experiments and modeling are well justified. However, in the case when population is a mixture of differently responding cells (for example apoptotic and proliferating or TNFα responding or not) the average trajectory does not represent any biological process and the model reproducing such trajectory is likely to be incorrect.
In Additional file 3 (Figs S1, S2, S3 and S4) we demonstrate the model's ability to reproduce major NF-κB pathway experiments on cells exhibiting oscillations in cytoplasmic to nuclear NF-κB localization that arise in response to persistent TNFα stimulation.
In the current paper we focus on TNFα signaling, a process initiated by binding of TNFα to the ubiquitous receptor TNFR1. In short, the action of regulatory pathway may be summarized as follows: Binding of TNFα trimer initiates receptor TNFR1 trimerization and formation of an active receptor complex in a multistep process involving binding of RIP and TRAF2. The active receptor complex activates the IKKK kinase (transformation from IKKKn to IKKKa). Active kinase IKKKa phosphorylates and activates the IKK kinase (transformation from IKKn to IKKa). Active IKKa kinase transiently binds to the cytoplasmic (NF-κB|IκBα) complex and phosphorylates IκBα initiating its degradation. Released NF-κB enters the nucleus to induce transcription of inhibitors IκBα and A20 genes. The first negative feedback loop involves the IκBα protein, which is rapidly resynthesized, enters the nucleus and recaptures NF-κB back into the cytoplasm. In the continued presence of IKKa, however, the resynthesized IκBα would be continuously degraded, which would result in continued nuclear NF-κB translocation. A second level of negative autoregulation occurs with the resynthesis of A20, a ubiquitin ligase which controls IKK activity. A20 initiates degradation of RIP, the key component of TNFR1 receptor complex, what attenuates the activity of receptors and directly associates itself with IKKa, enhancing its conversion to catalytically inactive IKKi. Inactive kinase IKKi spontaneously converts back to IKKn through the intermediate form IKKii. Similarly, active kinase IKKKa rapidly converts to the inactive form IKKKn.
Gene expression cascade (Fig. 2A) starts from the activation of a single gene copy, which may then serve as a template for the synthesis of tens or hundreds of mRNA molecules. In turn, a single mRNA molecule is a template for synthesis of hundreds of protein molecules. In this way the two IκBα gene copies are sufficient to replenish pool of IκBα proteins of about 100,000 molecules, within a half hour. As estimated experimentally by Yang  the endogenous level of IκBα molecules is 135,000, most of these molecules are degraded at in first 10 min. of high dose TNFα stimulation, and then IκBα level is approximately restored between 30 and 75 min. of stimulation .
At the level of cell surface receptors, a single TNFα trimer binding to the TNFR1 receptor leads to receptor trimerization and formation of a stable active receptor complex (Fig 2B). Grell  found that TNFα trimers dissociate from TNFR1 receptors with a half time of 33 min., while the internalization time is of order of 10 to 20 min. During this time the single active receptor may activate numerous IKKK kinase molecules. In turn, each active IKKKa activates numerous IKK kinases, and each of IKKa may phosphorylate several IκBα molecules leading to their degradation. The IKKK-IKK transduction cascade resembles the MAPK cascade and provides signal amplifica-tion of several orders of magnitude . This amplification mechanisms enables cells to respond to femtomolar concentrations of TNFα [17, 43–45] and references therein.
Recently, in cell population studies, Cheong et al.  observed activation of the NF-κB regulatory module in response to TNFα concentrations of 0.01 ng/ml. This equals to 200 fM (assuming that TNFα consists of trimers of mass 51 kDa) and implies about 1 TNFα trimer per 8 × 10-12 l or less than one TNFα trimer molecule per volume of mammalian cell which is of the order 2 × 10-12 l. This estimations suggests that cells may be activated by binding of a single, or few, TNFα trimers, and that at femtomolar concentration some cells may become active and some not, since TNFα binding is a stochastic process. It also indicates that when such small concentrations are considered, the average number of TNFα trimers per cell may be a better parameter to describe the experiment than the TNFα concentration itself. For example Chan and Aggarwal  observed two fold NF-κB induction at TNFα dose of 100 fM. For the EMSA assay they incubated 2 × 106 cells in 500 μl medium, which gives 5 TNF trimers per cell. The same concentration in tissue, where the cells are tightly packed would imply less than 1 molecule of activator per cell. The number of TNFR1 receptors per cell may vary significantly between cell lines , e.g. there are about 3000 TNFR1 per cell for HeLa , and about 10000 for Histiocytic lymphoma (U-937) cells , but much less for B-cell lymphoma (Raji) cells. Since in low dose experiments there are more TNFR1 receptors, than TNFα molecules, the concentration of free TNFα will be influenced by reaction with these receptors and will not remain constant in the course of low dose experiments. Similarly, when the spread of TNFα in the tissue is considered one may expect that TNFα diffusion will be strongly influenced by binding to free TNFR1 receptors, which may restrict cell to cell signaling to very short distances.
In HL60 cells, NF-κB activity was already observed at TNFα concentrations as low as 0.1 pM whereas maximum NF-κB activation required 0.4 to 2 pM TNFα, .
There is also an ample evidence that cells are able to respond to single viruses, which are known to activate NF-κB pathway through Toll-like receptors dependent and independent pathways, both engaging IKK, [46, 47]. Specifically, in a recent analysis of human A549 pulmonary type II epithelial cells infected by respiratory syncytial virus (RSV) at MOI = 1 (multiplicity of infection) we showed that 60% of cell exhibit RelA activation . The MOI = 1 implies (if the Poisson distribution of virions per cell is assumed) that 37% of cell will remain uninfected, while only 26% of cell will be infected by more than 1 virion. Thus the observed 60% fraction of responding cells implies that a single virus is enough to induce NF-κB activity in a cell. Arnold et al.  analyzed responses of peripheral blood mononuclear cells (PBMC) after exposure to low infectious RSV doses, with MOI from 0.001 to 1. Even at MOI as low as 0.01–0.1 they observed pronounced secretion of, NF-κB responsive cytokines like IL-8, IL-6 and TNFα at 4 hours after infection.
The experiment helped us to determine the parameters governing the NF-κB – A20 – IKK negative feedback loop, Figs 3G and 3H. Although after a short break the system did not respond, full recovery of the system was observed after 2.5 hours. Surprisingly, when the duration of the break was extended to 3 hours the second peak of IKK activity was higher than the first one. One could interpret this result as the evidence that after 2.5 to 3 hours the cytoplasmic A20 concentration decreases to pretreatment values and some other protein is elevating IKK activity. One of the candidates could be TRAF2, which is NF-κB responsive and is a constituent of the TNFR1 receptor complex . To test this hypothesis we assume that at second peak the activity of TNFα bound receptors is elevated 10 fold, Fig. 3I, and in fact this modification resulted in an elevated IKK profile similar to the experimental findings. However, since we do not find this a strong enough verification, we have not introduced this modification to the current model.
Mentioned earlier, the ability of the cell to respond to femtomolar TNFα concentrations suggests that NF-κB regulatory module may be turned on by the single TNFα trimer activating one TNFR1 receptor (receptor trimer). According to our model the lower TNFα dose results in lower probability of cell activation, but not in a much weaker response in the responding cells (Figs. 5A and 6B). We have chosen the receptors activation coefficient such that 90% of cells are activated (have at least one receptor active) in the first 10 minutes of TNFα stimulation by the dose of 1 ng/ml, which may be treated as a saturation dose, until the persistent stimulation is considered. This is in agreement with experimental observations in which responses to TNFα doses of 1 ng/ml or 10 ng/ml are almost identical. Below this dose, the individual cell responses become asynchronous and at any given moment the population is a mixture of responsive and nonresponsive cells.
Because the NF-κB regulatory network has the amplification-saturation pathway build in, the final cell response measured at the level of NF-κB responsive genes is about the same in the case when the single receptor (trimer) is activated as when 100 receptor complexes are activated. To see how cell's response depends on the number of active receptors, one can perform the following numerical experiment. Let us assume that at t = 0, a given number N (1, 3, 10, 30 or 100) of receptors is activated and that these receptors remain active for 10 min. We traced the average heights of the peaks of active IKKK, active IKK, nuclear NF-κB, and A20 mRNA, Fig. 5B. We thus observed four levels of signal saturation; first three are connected with the transduction-amplification pathway, the last one with gene regulation. According to our model the fivefold decrease in peak IKK activity when we pass from 100 to 1 active receptor, does not result in much lower NF-κB oscillation amplitude. This is in agreement with bulk of experimental data suggesting that peak IKK activity at high and intermediate TNFα dose of 10 or 1 ng/ml is much higher than needed for efficient IκBα degradation. Specifically, Lee et al.  demonstrated (for TNFα dose of 10 ng/ml) that a very low IKK activity tail is sufficient to maintain oscillations in IκBα level, and that a higher IKK activity tail, typical for A20 deficient MEFs, totally suppress IκBα accumulation. This result was confirmed by quantitative data of Werner et al.  who stimulated A20 deficient, immortalized 3T3 cells by a 45 minute long TNFα pulse at a dose of 1 ng/ml. They found that the peak of IKK activity (same for wild type and A20 -/- cells) was followed by tail of fivefold lower magnitude but that this was enough to keep most of NF-κB in the nucleus. The fact that NF-κB remains nuclear implies high IκBα synthesis rate and thus high degradation rate, so it may not accumulate to uptake NF-κB back to cytoplasm.
The high sensitivity of NF-κB responses to low dose TNFα stimulation has been already reported by Cheong et al.  who, based on experiments with decreasing TNFα dose from 10 to 0.01 ng/ml, concluded that NF-κB amplitude is a logarithmic function of TNFα dose: It decreases only twice as the doses changes from 10 to 0.1 ng/ml. The effect we reported in Figs. 5B and 6B is even more dramatic, the NF-κB amplitude remains almost independent of the TNFα dose assuming that only the activated cells are considered. The single cell model we constructed suggests that the decrease of NF-κB amplitude shown by Cheong et al.  (Fig. 1) is at least partially due to pure synchronization of cells. As recently experimentally demonstrated by Sillitoe et al.  the averaging causes that NF-κB pulses appear much smaller than they are. Sillitoe et al.  experiment was performed for TNFα dose of 10 ng/ml, but even at that dose (at which one may expect cells are relatively well synchronized) the first peak, amplitude of the average trajectory is about twice smaller than single cell amplitudes, and this difference is larger for the second peak at which cells are less synchronized. In addition, at the smallest dose of 0.01 ng/ml studied by Cheong et al. , which corresponds to only 1/4 of TNFα trimer per cell volume, we observe the effect of averaging over responsive and nonresponsive cells, which lowers the amplitude of average NF-κB trajectory.
There is still not enough experimental data at single cell level to uniquely determine how each step of transduction amplification cascade contributes to amplification. Cheong et al.  attempted to address this question based on population data and found that a 100 fold decrease of TNFα dose from 10 to 0.1 ng/ml results in a fourfold decrease of IKK activity peak and in a twofold decrease of NF-κB amplitude. As already discussed, the population data can be quite misleading, and in fact to fit to this data Cheong et al.  had to assume time dependent IKK activation and inactivation rates. Such approach seems artificial and cannot be applied to more general protocols of TNFα stimulation (e.g., to pulse-pulse stimulation).
As discussed above, the constructed model has the property that the stable activity of a single receptor is sufficient to turn on the NF-κB regulatory module and trigger the transcription of NF-κB dependent genes. Although the experimental data suggests that the a very limited number of active receptors is needed to induce NF-κB activity, there is no evidence that activity of a single receptor is sufficient. One may expect that there exists a threshold of the number of active receptors (receptor complexes) needed for cell activation. In the case of T-cells stimulation, it was experimentally demonstrated  that just three major peptide histocompatibility complexes are required to induce T-cell cytotoxic activity. T-cell receptor signalling involves MAPK kinase cascade , which can convert graded inputs into switch-like outputs . Both MAPK and MAPKK require two phosphorylations to become fully activated, and this dual phosphorylation is a source of nonlinearity in signal processing, which together with saturation (with a strong signal all MAPK is phosphorylated) results in a switch-like output.
Since we did not identify any simple mechanisms defining the threshold here, we followed the simplest possibility and assumed that amplification provided by the TNFR1-IKKK-IKK-IκBα cascade is high enough to cause that activation of a single receptor is sufficient to induce cell activity. To rule out or confirm the competitive threshold hypothesis it is worth to consider a series of single cell experiments, in which cells would be stimulated for 5 min. with a decreasing TNFα dose, and measure the fraction of responding cells. Pulse stimulation would help localize activation and response to relatively short time periods.
According to our model, the probability that a single receptor is activated by a TNFα stimulation lasting for time t at a concentration C is equal to
pr1(C) = 1 - exp(-k b Ct), (1)
The single pulse, low dose, experiment would also enable us to determine the "minimum cell response" as the average (over subpopulation of responding cells) NF-κB oscillation amplitude for very small C. According to our model (Fig. 6B) such a minimum response exists and moreover it is high enough to assure unequivocal expression of NF-κB responsive genes in majority of responding cells.
Stochastic gene activation (leading to the burst of proteins) and stochastic cell activation (leading to the massive NF-κB nuclear translocation) leads to what might be called "stochastic robustness" in cell regulation. If a given gene is activated, a large burst of proteins is produced, in order to assure a sufficient level of activity of these proteins. Stochastic robustness assures the minimal response to the signal. Decreasing magnitude of the signal mostly reduces the probability of response, which leads to a smaller fraction of responding cells. This may be a useful strategy: If the TNFα signal is low, some cells respond by a massive NF-κB translocation, whereas some do not respond at all. It helps to avoid ambiguity, such as when a small nuclear concentration of NF-κB leads to activation of an undefined fraction of NF-κB responsive genes. Such an ambiguous response might do more harm than good. In fact the all-or-none responses arising due to ultrasensitivity and saturation or bistability (typically connected with positive feedbacks) have been reported for various signalling elements: TCR signaling , Xenopus p42 MAPK cascade , transduction cascade JNK in oocytes , and Lac operon, e.g. . The growing evidence of bistability in various system may suggest that it is a good strategy for regulation at the tissue level.
Stochastic robustness allows cells to respond differently to the same stimulation, but makes their individual responses better defined. Both effects could be crucial in early immune response: Diversity in cell responses causes the tissue defense to be harder to overcome by relatively simple programs coded in viruses and other pathogens. The more focused single cell responses help cells to decide their individual fates such as proliferation or apoptosis.
The proposed model supports hypothesis that at nanomolar concentrations NF-κB activity results from binding of single TNFα trimers. Such hypersensitivity of NF-κB regulatory module may be the only way to detect and respond to single viruses invading cells and to allow for cytokine extracellular signaling.
The model applied involves two-compartment kinetics of NF-κB and its inhibitors A20 and IκBα, and allows analyzing single cell responses to TNFα stimulation of arbitrary and time dependent intensity, in both wild type and A20 deficient cells. It combines our previous model  with the signal transduction-amplification cascade, which transmits the signal from the TNFR1 receptors, , to the NF-κB-IκBα-A20 regulatory module. The single event of receptor activation is amplified by the three-step signal transduction cascade involving activation of kinase named IKKK (IKK activating kinase) and IKK . Biologically, there are at least two kinases involved in this process: MEKK3 [58, 59] and TAK1 [59–61]. Shim et al.  showed that TAK1 mediates IKK activation in TNFα and IL-1 signaling pathways. Homozygous mutants (Tak1m/m) do not show NF-κB activity under TNFα stimulation, and exhibit lowered NF-κB activity than the wild type cells in response to IL-1. Nho et al.  found that silencing MEKK3 by siRNA reduces IKK activity. These two findings and the experiment by Blonska et al.  suggest that TAK1 is recruited to the TNFR1 complex via RIP and likely cooperates with MEKK3 to activate NF-κB in TNFα signaling. In our model we mimic this possibly complex transduction mechanism by a single entity: IKKK. To accomplish the above, we assume that IKKK migrates toward the receptor and is activated at a receptor (transformation from IKKKn to IKKKa). Active IKKKa molecules activate IKK molecules (transformation from IKKn to IKKa) which in turn phosphorylate IκBα molecules leading to their ubiquitination and degradation. It is assumed that the total number of IKKK and IKK molecules (as well as of the NF-κB molecules) is constant; i.e., their degradation is balanced by production, but both terms are omitted in the mathematical representation. The IKKK molecules may exist in one of two states,
native, neutral, IKKKn, specific to unstimulated resting cells without, and
The activity of IKKK is very transient, i.e., it is activated rapidly (with the maximum rate 1/s from single active receptor) and inactivated with the half time of about 1 min.
IKK complexes, consisting of catalytic subunits IKKα and IKKβ and regulatory subunits IKKγ, may exist in one of four states,
native neutral (denoted by IKKn), specific to unstimulated resting cells.
active (denoted by IKKa), arising from IKKn via phosphorylation of serines 177 and 181 of IKKβ subunits  upon the IKKKa induced activation,
inactive (denoted by IKKi), but different from the native neutral form, arising from IKKa possibly due to overphosphorylation, and
transient between IKKi and IKKn, also inactive, denoted by IKKii.
The activation pathway considered, which involves kinases IKKK and IKK, resembles the one introduced by Park et al. . The main difference is that we take into account the inactivation of IKKa is due to its transformation to the state IKKi different from the native state IKKn. As a result, in contrast to Park et al. model , the tonic TNFα stimulation induces transient IKK activity in agreement with the experimental data, Figs 6 and S2.
In resting cells, the unphosphorylated IκBα binds to NF-κB and sequesters it in an inactive form in the cytoplasm. IKKa mediated phosphorylation of IκBα leads to it degradation and releases the main activator NF-κB, which then enters the nucleus and triggers transcription of the two inhibitors and numerous other genes. The newly synthesized IκBα leads NF-κB out of the nucleus and sequesters it in the cytoplasm.
IKK inactivation is controlled by the second inhibitor A20, which like IκBα; is strongly NF-κB responsive . The exact mechanism of A20's action is still not fully resolved. Here we assume that A20 acts in two ways: (1) It initiates degradation of RIP, the key component of the TNFR1 receptor complex , which attenuates the activity of receptors, and (2) it directly associates with IKKγ , enhancing IKKa conversion (this conversion takes place also in A20 deficient cells but at a slower rate) to catalytically inactive form IKKi. The exact mechanism of IKK inactivation remains unresolved: According to Delhase et al.  IKK inactivates via autophosphorylation of serines in IKKβ C-terminal region. However, recently Schomer-Miller et al.  found that this autophosphorylation does not diminish IKK activity and suggested that phosphorylation of serines 740 and 750 in NBD/γBD domain of IKKβ may have a regulatory role and that their phosphorylation may downregulate IKK activity. The form IKKi spontaneously converts into IKKn through inactive intermediate forms collectively denoted by IKKii. The number of these forms may be large since there are at least 16 serine residues in IKKβ , which may be involved in regulation of IKK activity. This intermediate step is introduced in the model to account for the delay needed to process the inactivated form IKKi into native state IKKn. This effect is manifested in A20-/- cells (at persistent or long lasting TNFα stimulation) as a downregulation of IKK activity at about 30 min. followed by a higher plateau (Figs. S2 and S3, reproduced after  and  experiments). According to our model in the first minutes of high dose TNFα stimulation most of IKKn is used up so the IKK activation rate is low, only after some IKKn is recovered via intermediated form IKKii, the activation rate and the level of IKKa may increase.
The inhibitor IκBα migrates between the nucleus and cytoplasm and forms complexes with NF-κB molecules. The nuclear IκBα|NF-κB complexes quickly migrate into the cytoplasm. The second inhibitory protein A20 is considered only in the cytoplasm where it triggers the inactivation of IKK. It is assumed that the transformation rate from IKKa into IKKi is the sum of the constant term and a term proportional to the amount of A20.
The transcriptional regulation of A20 and IκBα genes is governed by the same rapid elongation regulatory mechanism with a rapid coupling between NF-κB binding and transcription. The mechanisms for NF-κB dependent regulation of IκBα and A20 are based on the control of transcriptional elongation. In this situation, stalled RNA polymerase II is rapidly activated by NF-κB binding to enter a functional elongation mode, and requires continued NF-κB binding for reinitiation. This is represented in our model by tight coupling of NF-κB binding to mRNA transcription. We assume that all cells are diploid, and both A20 and IκBα genes have two potentially active homologous copies, each of which is independently activated due to binding of NF-κB molecule to a specific regulatory site in gene promoter. Following [1, 3, 8, 35] and others we made the simplifying assumption that each gene copy may exist only in one of two states; active and inactive. When the copy is active the transcription is initiated at a high rate, when the copy is inactive transcription is inhibited. The gene copy becomes inactive when the NF-κB molecule is removed from its regulatory site due to the action of IκBα molecules, which bind to DNA-associated NF-κB, exporting it out of the nucleus.
In this work, as in our recent papers  we follow the method proposed by Haseltine and Rawlings  and split the reaction channels into fast and slow. We consider all reactions involving mRNA and protein molecules as fast and the reactions of receptors and genes activation and inactivation as slow. Fast reactions are approximated by the deterministic reaction-rate equations, whereas slow reactions are considered stochastic.
According to the above, the mathematical model consists of 15 ordinary differential equations (ODEs) accounting for
formation of the (IκBα|NF-κB) complexes,
IKKK and IKK kinase activation and inactivation
IKKa driven IκBα phosphorylation,
A20, IκBα and phospho-IκBα proteins degradation
transport between nucleus and cytoplasm, and
transcription and translation.
All the substrates are quantified by the numbers of molecules. The upper-case letters denote substrates or their complexes. Nuclear amount is represented by subscript n, while subscript c denoting amount of substrate in the cytoplasm is omitted, to simplify the notation. Amounts of the mRNA transcript of A20 and IκBα are denoted by subscript t:
IKKn – neutral form of IKK kinase,
IKKa – active form of IKK,
IKKi – inactive form of IKK,
IKKii – inactive intermediate form of IKK,
K NN – total number of IKK = IKKn+IKKa+IKKi+IKKii molecules (assumed to be constant in time)
IKKKa – amount of active form of IKKK,
IKKKn – amount of neutral form of IKKK,
K N – total number of IKKK=IKKKn+IKKKa molecules (assumed to be constant in time)
IκB – cytoplasmic amount of IκBα,
IκB n – nuclear IκBα,
IκB t – IκBα mRNA transcript,
IκB p – phosphorylated cytoplasmic IκB
NFκB|IκB – cytoplasmic (NF-κB|IκBα) complexes
NFκB|IκB p – phosphorylated cytoplasmic IκBα complexed to NFκB
NFκB|IκB n – nuclear (NFκB|IκBα) complexes
TNF – TNFα concentration,
G I κ B – discrete random variable, state of IκBα gene,
GA20 – discrete random variable, state of A20 gene,
k v = V/U – ratio of cytoplasmic to nuclear volume.
B – number of active receptors, M – total number of receptors (assumed to be constant in time)
Note that Eq. (15) (as well as 12) naturally produces saturation in transcription speed. When the nuclear amount of regulatory factor NF-κB is very large, then the binding probability is much larger than the dissociation probability, and the gene state would be G a = 2 for most of the time. In such case the transcription would proceed at a maximum rate, 2c 1.
We assume that both A20 and IκBα genes have two homologous copies independently activated due to NF-κB binding, and inactivated due IκBα mediated removal of NF-κB molecules, and that binding and dissociation propensities r b (t) and r d (t), respectively, are equal for each copy:
r b (t) = q1 × NFκB n (t), r d (t) = q2 × IκBα n (t).(19)
The state of gene copy G i (i = 1, 2) is G i = 1 whenever NF-κB is bound to the promoter regulatory site, and G i = 0 when the site is unoccupied. As a result the gene state G = G 1 + G 2 can be equal to 0, 1 or 2. In this approximation the stochasticity of single cell kinetics solely results from discrete regulation of receptors activity and transcription of A20 and IκBα genes.
In model computations, the amounts of all the substrates are expressed as the numbers of molecules. Since we use the ODE's to describe most of the model kinetics, amounts of molecules are not integer numbers, but since these numbers are in most cases much greater than 1, such description is reasonable.
The numerical scheme implemented follows that of :
(2) We select two random numbers p 1 and p 2 from the uniform distribution on [0, 1].
(4) There are 6 potentially possible different reactions:
receptor may be activated or inactivated. Typically in time course there are many inactive receptors which may activated and active receptors which may be inactivated, but since the receptors are assumed to be identical it is not important which one of them changes its state.
NF-κB may bind to or dissociate from any of two alleles of A20 and IκBα genes.
where r i (t + τ), i = 1,...,6 are individual reaction propensities and k is the index of the reaction to occur.
(5) Finally time t + τ is replaced by t, and we go back to item (1).
In all simulations before TNFα stimulation starts, we simulate a resting cell for time t randomly chosen from the interval of 15 to 25 hours in order to get equilibrated and randomized initial conditions. As shown in  the resting cells oscillate. Due to natural degradation of IκBα, some NF-κB molecules may occasionally enter the nucleus and activate the A20 or IκBα gene, which results in bursts of A20 and IκBα mRNAs and proteins. As a result the initial (at t = 0) level of A20 or IκBα mRNA and protein differs across the population, which influences future cells evolution.
To compare model predictions with Nelson et al. experiment  we calculate the NF-κB nuclear to cytoplasmic oscillation amplitudes assuming that 40% of total pool of NF-κB does not oscillate, Fig. 4. We also convert ratio of nuclear to cytoplasmic amounts to ratio of light intensities, making use of the fact that nucleus is smaller than cytoplasm. These two operations do not influence the character of the oscillations but only their amplitude. The non-oscillatory fraction (not taken in account in the model) remains in the cytoplasm even after all IκBα is degraded. This fraction is observed also in population experiments (data not shown) and it possibly consists of NF-κB bound to other IκBα isoforms and of NF-κB molecules with nuclear localization sequence damaged.
The model is intended. cells showing oscillations in NF-κB nuclear-to-cytoplasmic localization under persistent TNFα stimulation. We do not to fit the parameters to any given single experiment but we intuitively choose the set of parameters which produces qualitative agreement with a major subset of existing data. Quantitative agreement is not possible because there is a substantial discrepancy mostly in IKK activity and NF-κB peak timing between experiments). Because of a large number of undetermined parameters, this is a tedious task, but in our opinion it is better to produce a model in qualitative agreement with the current and previous experiments, than a model perfectly fitted to a single experiment with limited set of data. We based our choice of parameters on both single cell and population experiments. In the letter case, to compare our model with experiments, we average a large number of single-cell stochastic simulations. This procedure is much more time consuming that comparing the deterministic model with the population data, but as already shown (in the case of low dose TNFα stimulation) the population data do not correspond to any biological process, and thus constructing the model fitted to such a data is not justified.
We applied the following method of choosing the values of parameters:
Start from a reasonable set of parameters, which produces a correct steady state in the absence of TNFα signal.
Proceed with the signal initiated by TNF downstream the autoregulatory loops.
Iterate item 2 until the fit to all the data is satisfactory.
The TNFα signal first causes activation of receptors then transformation of IKKKn into IKKKa which catalyses transformation of IKKn into IKKa. In turn, IKKa catalyses degradation of cytoplasmic (IκBα|NF-κB); enabling the free NF-κB to enter the nucleus. Once NF-κB builds up in the nucleus it upregulates the transcription of the A20 and IκBα genes. After being translated, A20 facilitates transformation of IKKa into IKKi and blocks the receptors, while IκBα enters the nucleus, binds to NF-κB and leads it into cytoplasm. As stated in item 2, we first fit the coefficients regulating IKK activation (using data on IKK activity), then the coefficients regulating degradation of the cytoplasmic (IκBα|NF-κB) and IκBα degradation and so forth. If there were no feedback loops in the pathway, the proposed method would be quite efficient, but, since they exist, it is necessary to iterate the signal tracing several times, until the fit is satisfactory. Once a satisfactory set of parameters is found, we observe that this set of parameters is not unique. This ambiguity is mainly caused by the lack of measurements of absolute values of protein or mRNA amounts. The action exerted by some components of the pathway onto the rest of the pathway is determined by their amounts multiplied by undetermined coupling coefficients. Hence, once we have a good set of parameters, we may obtain another one using a smaller coupling coefficient and by proportionately enlarging the absolute level of the component. Since not all parameters may be determined based on existing data we have assumed values of part of parameters mostly based on our intuition and fitted the remaining ones. By such approach we show how much information can be inferred from available experimental data. Ambiguity in parameter determination leads to significant differences between parameters of our model and the corresponding parameters chosen by others groups of researcher. Values of all model parameters are listed in discussed in detail in Additional file 2.
The validation of the proposed model, is based on our data (see Fig. 3) in addition to  (see Fig. 4),  (see Fig. 7),  (see Additional file 3: Fig. S1),  (see Additional file 3: Fig. S2),  (see Additional file 3: Fig S3),  (see Additional file 3: Fig S4) and  experiments. Typical experimental data at our disposal consist of measurements made at time points that are not uniformly distributed. Therefore, to compare our predictions to experimental data, we rescale the time coordinate, i.e. we calculate the values of corresponding functions in the same time points as in the experiment, and then to guide the eye we connect the resulting discrete points by straight-line segments, thus obtaining a saw-like graph.
TNF alpha inducible protein 3
inhibitor of κB α subunit
interleukin-1 receptor-associated kinase 1
mitogen-activated protein kinase kinase kinase 3
receptor interacting protein
transforming growth factor-β-activated kinase 1
tumor necrosis factor α
tumor necrosis factor receptor associated factor 2.
We thank MRH White and members of his laboratory for discussion and critical reading of the manuscript. The MATLAB program written to perform the simulations will be available at our website .
This work was supported by Polish Committee for Scientific Research Grants No. 4T07A 001 30, 3T11A 019 29 and PBZ-MNiI-2/1/2005 and by NHLBI contract N01-HV-28184, Proteomic technologies in airway in inflammation (A. Kurosky, P.I.).
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.