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.
Results and discussion
Stochastic switches and amplification cascades in NF-κB regulation
Our considerations are based on the two-feedback loop stochastic model of the NF-κB pathway, which combines the signal transduction cascade that connects cell surface receptors with the core regulatory module analyzed previously . Current model involves two-compartment kinetics of transcription factor NF-κB, its activators, IKKK, IKK and inhibitors, A20 and IκBα, shown in Fig. 1. IKKK represents the IKK activating kinase, which itself is activated at the TNFR1 receptor complex (see Materials and Methods for details).
Figure 1. Two feedback-loop model of NF-κB regulatory pathway.
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.
Additional file 3. Supplementary figures S1, S2, S3 and S4. Validation of the proposed model based on experiments by Hoffmann et al. , Fig. S1; Lee et al. 2000 , Fig. S2; Werner et al. , Fig S3; and Nelson et al. , Fig S4.
Format: PDF Size: 2.5MB Download file
This file can be viewed with: Adobe Acrobat Reader
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.
We identified two stochastic processes crucial to the functioning of the NF-κB regulatory pathway: (1) Activation of A20 and IκBα genes via binding of NF-κB molecules to the genes promoters and (2) activation of TNFR1 receptors via binding of TNFα trimers. These stochastic events may influence the evolution and fate of the cell due to their association with amplification cascades, as shown in Fig. 2 (see [38-40] for discussion of signal and noise propagation in gene neworks and amplification cascades).
Figure 2. Two ways of signal amplification present in the model. Panel A: A20 (or IκBα) gene expression. Panel B: Transduction amplification pathway leading to IκBα degradation.
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-12l or less than one TNFα trimer molecule per volume of mammalian cell which is of the order 2 × 10-12l. 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.
System recovery: pulse-pulse experiment
Although several studies [23,25,45,50,51] have clarified the relationship between time profile of IKK activity and NF-κB oscillations, less is known about the regulation of IKK activity itself. In response to TNF stimulation, IKK is rapidly activated and then inactivated. As demonstrated by the experimental findings  and modeling , an important mechanism for this inactivation is mediated by the NF-κB – A20 negative feedback loop. Here, a short pulse of TNFα results in a brief peak of IKKa followed by peak in nuclear NF-κB and burst of A20 mRNA and protein. The newly synthesized A20 attenuates IKK activity, making the cell temporarily less sensitive to subsequent TNFα stimulation. As described earlier, A20 acts via direct and indirect mechanisms. Directly it binds to the regulatory subunit of IKK, speeding transformation of IKKa into the inactive, hyperphosphorylated IKK (represented as IKKi). Indirectly, it ubiquitinates RIP, a necessary component of the active receptor complex leading to its specific proteasomal degradation, thereby lowering the average receptor activity. An interesting question arises: How much time is needed for the system to recover, i.e., how long should be the break after a brief TNFα pulse to allow for the equal response to subsequent pulse? In order to answer the question we performed an experiment (see Additional file 1 for details of experimental protocol) in which the population of 3T3 cells is stimulated by two brief TNFα pulses (5 min., at 20 ng/ml) spaced at various times apart ranging from 30 to 180 min. This saturating TNFα concentration makes the individual cell responses well synchronized, allowing for reliable population analysis, Fig 3. The IKK kinase activity is measured 5 min. after the second pulse, and the IKK activity of the second peak relative to the first is expressed as a function of break duration.
Format: PDF Size: 26KB Download file
This file can be viewed with: Adobe Acrobat Reader
Figure 3. Pulse-Pulse experiment. Panels A to F: IKKa and NF-κB oscillation amplitude resulting from single cell numerical simulations. Two TNFα (20 ng/ml), 5 min pulses are separated by 30', 60', 90' 120' 150' or 180' breaks. Four single cell simulations (marked by different colors) are shown in each panel. Panel G: IKK kinase activity in 3T3 cells stimulated by two 5 min. pulses of TNFα (20 ng/ml), separated by 30', 60', 90' 120' 150' or 180' break. First sub-blot corresponds to unstimulated cells, second sub-blot shows intensity of first peak, while next six sub-blots show the intensity of second peak as a function of brake duration, see Additional file 1 for the protocol of the experiment. Panels H and I magnitude of the second and the first peak of IKK and NF-κB activity (nuclear NF-κB and IKKa) calculated from average over 500 cells. In numerical experiment presented in Panel I, the receptor activity coefficient ka was elevated 10 fold before second TNFα pulse.
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.
Cell respond to stimulations by a wide range of TNFα doses
We performed the single cell stochastic numerical simulations (SCSNS) of our model to analyze the individual cell responses to persistent stimulation in a broad range of TNFα doses. We assumed that there were 1000 TNFR1 receptors on cell surface, and each of these receptors might be independently activated or inactivated, and that the activation rate was proportional to the TNFα concentration, which was kept constant during the experiment. At a high TNFα dose (above 1 ng/ml) the receptor activation rate is high and most of cells are activated in the first few minutes after the TNFα stimulation begins. As a result, the first peak of IKK activity and the first peak of NF-κB nuclear to cytoplasmic oscillation are well synchronized among cells. Synchronization of subsequent peaks of NF-κB oscillations decreases due to the stochastic processes of activation of TNFR1 receptors and A20 and IκBα genes. Such behavior well agrees with Nelson et al.  single cell experiment on SKNAS cells, Fig 4A versus 4B.
Figure 4. Single cell experiment and simulation. Panel A: Nuclear-to-cytoplasmic NF-κB oscillations measured by Nelson et al.  shown for four SKNAS cells (marked by different colors) stimulated by 10 ng/ml TNFα. Panels B to F: IKKa and NF-κB oscillation amplitude resulting from model simulations, respectively for TNFα dose: 10, 1, 0.3, 0.1 and 0.03 ng/ml. In each panel we show three single cell simulations, marked by different colors.
For low doses (0.1 and 0.03 ng/ml, Figs 4E and 4F), the activation of each cell is typically due to the activation of single receptor or a small number of receptors and thus the first response time varies between cells. As a result, the NF-κB oscillations are not synchronized at all. Moreover, for the lowest dose (0.03 ng/ml, Fig 4F) the time span between subsequent oscillations varies significantly. The last is due to the fact that when TNFα molecules dissociate the cell become unstimulated. It is however possible that dissociated TNFα molecule is immediately captured by other TNFR1 receptors (or transiently by TNFR2 receptors) of the same cell. This effect would keep the once stimulated cell in semi-periodic oscillatory mode for the longer time. Interestingly, the amplitude of the second and subsequent oscillations is higher for the low TNFα dose than for the high dose. This is due to the fact that at a low TNFα dose the IKK activation proceeds not as fast: the first peak of IKKa is smaller than for a high dose (Fig. 5) and more IKKn remains to be activated. As a result the system at low dose stimulation exhibits subsequent IKKa pulses followed by high NF-κB oscillations (Fig. 4).
Figure 5. Characteristics of cell activation as a function of dose. Panel A: Probability of cell activation (at least one receptor activation) for 10 min. stimulation as a function of TNFα dose. Panel B: Cell activity measured as an average height of the IKKKa, IKKa, nuclear NF-κB, and A20mRNA peak resulting from 10 min. activity of 1, 3, 10, 30, or 100 receptors.
Activation of a single TNFR1 receptor complex (trimer) turns on the NF-κB regulatory module
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
where kb is the receptor activation rate. Thus, if activation of a single receptor is sufficient to activate the cell, then the fraction of responding cells F would decrease approximately linearly with TNFα concentration C for low concentrations C. If however, at least n receptors must be activated in order to activate the cell then the probability of cell activation would be equal to prn
where N is total number of receptors and kbn is the new hypothetical TNFα binding constant. In such case the fraction of responding cells F would decrease faster, approximately as Cn (for small C). In Fig. 6A we compare model predictions (n = 1, red line) against threshold hypothesis (n = 3, green line) and (n = 10, black line) for N = 1000. Activation probability functions pr3(C) and pr10(C) are normalized by adjusting activation rates kb3 ≃ 5.03 kb and kb10 ≃ 20.79 kb, what assures that all prn(C) functions intersect at the same point at which prn(C) = 1/3.
Figure 6. Model predictions for pulse low dose TNFα stimulation. Panel A: Probability of cell activation during 5 min. long TNFα stimulation as a function of dose, model prediction (red) versus threshold hypothesis: green (3 receptors), black (10 receptors). Panel B: Average height of nuclear NF-κB peak calculated for subpopulation of N = 500 of activated cells (at least with one receptor activated) after 5 min. long TNFα stimulation as a function of TNFα dose.
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 stochastic model proposed by us explains the mechanisms of TNFα cell activation at nanomolar concentrations when the number of molecules of activation factor per cell is very limited. In the example of TNFα diffusion in the tissue, considered by Cheong et al. , concentration of 0.01 ng/ml or 200 fM corresponds to less than one TNFα trimer per cell. In such case population consists of responding and nonresponsive cells, and the output observed at the population level reflects averaging these two subpopulations. At low dose stimulation no individual cell evolves like the average and the average does not correspond to any biological process. This means that the model build to follow the average trajectory in general may not be correct. Despite this fact, the experimental average may be useful when compared with the average of many single cell stochastic simulations. Predictions from our model agree qualitatively with large set of population data (see Additional file 3: Figs S1, S2, S3 and S4) including the Cheong et al.  experiment in which IKK and NF-κB activity were measured across a wide range of TNFα concentrations, Fig. 7.
Figure 7. Average IKK and NF-κB activity (nuclear NF-κB and IKKa) as function of TNFα dose, model versus Cheong et al. 2006 experiment . Panels A: IKKa and IKK kinase activity, Panel B: nuclear NF-κB calculated as average of N single cell simulation: N = 50 for TNFα dose equal 10 ng/ml, N = 100 (1 ng/ml), N = 500 (0.1 ng/ml) and N = 2500 (0.01 ng/ml). and EMSA assays of nuclear NF-κB. 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. Despite the qualitative agreement between the model and the experiment shown there is substantial difference in IKK activity peak timing and smaller difference in NF-κB profile timing for 10 ng/ml dose. In fact there is also large discrepancy between different experiments; here the IKK activity peak is observed for 5 min., while in  for 10 min. and in  for 15 min.
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.
Stochastic NF-κB model with two feedback loops
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
• active IKKKa.
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,
• KNN – 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,
• KN – total number of IKKK=IKKKn+IKKKa molecules (assumed to be constant in time)
• IκB – cytoplasmic amount of IκBα,
• IκBn – nuclear IκBα,
• IκBt – IκBα mRNA transcript,
• IκBp – phosphorylated cytoplasmic IκB
• NFκB|IκB – cytoplasmic (NF-κB|IκBα) complexes
• NFκB|IκBp – phosphorylated cytoplasmic IκBα complexed to NFκB
NFκB|IκBn – nuclear (NFκB|IκBα) complexes
• TNF – TNFα concentration,
• GIκB – discrete random variable, state of IκBα gene,
• GA20 – discrete random variable, state of A20 gene,
• kv = V/U – ratio of cytoplasmic to nuclear volume.
• B – number of active receptors, M – total number of receptors (assumed to be constant in time)
IKKK in active state (IKKKa)
The first term describes IKKK kinase activation, i.e. transformation from IKKKn (amount of which is IKKKn = KN - IKKKa) due to action of active receptors B(t), whose activity is attenuated by A20. The second term describes spontaneous inactivation of the kinase
IKK in the natural state IKKn
The first term describes IKKn recovery from the intermediate form IKKii (amount of which is IKKii = KNN - IKKn - IKKa - IKKi), whereas the second term describes depletion of IKKn due to its transformation into IKKa mediated by IKKKa
IKK in the active state IKKa
The first term represents transition from IKKi to IKKa mediated by IKKKa, whereas the second term represents depletion of IKKa due to its transformation into inactive form IKKi mediated by A20
IKK in the inactive state IKKi
The first term corresponds to the formation of inactive IKKi from IKKa by A20 mediated inactivation, whereas the second term describes transformation into IKKii
The first term describes IκBα phosphorylation due to catalytic action of IKKa, the second term catalytic degradation of phosphorylated IκBα
Phospho-IκBα complexed to NF-κB (NFκB|IκBp)
The first term describes IκBα phosphorylation (in complexes with NF-κB) due to the catalytic action of IKKa, the second term catalytic degradation of phosphorylated IκBα (NF-κB is recovered)
Free cytoplasmic NF-κB
The first two terms represents liberation of free NF-κB due to degradation of IκBα in (IκBα|NF-κB) complexes and its depletion due to formation of these complexes. The third term accounts for liberation of NF-κB due to degradation of phospho-IκBα. The last term describes transport of free cytoplasmic NF-κB to the nucleus,
Free nuclear NF-κB
The first term describes transport into the nucleus. The second term represents depletion of free nuclear NF-κB due to the association with nuclear IκBα and is adjusted, by multiplying the synthesis coefficient a1 by kv = V/U, to the smaller nuclear volume resulting in a larger concentration,
Described by its mRNA synthesis and constitutive degradation,
The first term stands for NF-κB inducible synthesis, while the second term describes degradation of the A20 transcript,
Free cytoplasmic IκBα protein
The first term accounts for IKKa induced phosphorylation, the second for NF-κB binding. The second line describes IκBα synthesis and the constitutive degradation of IκBα. The last two terms represent transport into and out of the nucleus,
Free nuclear IκBα protein
The first term corresponds to IκBα association with nuclear NF-κB (adjusted, by multiplying the synthesis coefficient a1 by kv, for the smaller nuclear volume resulting in larger concentration), and the last two terms represent the transport into and out of the nucleus,
The first term stands for NF-κB inducible synthesis, whereas the second term describes degradation of IκBα transcript
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 Ga = 2 for most of the time. In such case the transcription would proceed at a maximum rate, 2c1.
Cytoplasmic (IκBα|NF-κB) complexes
The first line describes formation of the complexes due to IκBα and NF-κB association and their degradation. The first term in the second line represents phosphorylation of the (IκBα|NF-κB) complexes due to the catalytic activity of IKKa. The last term represents transport of the complex from the nucleus,
Nuclear (IκBα|NF-κB) complexes
Described by their formation due to IκBα and NF-κB association (adjusted, by multiplying the synthesis coefficient a1 by kv, to the smaller nuclear volume resulting in larger concentration) and their transport out of the nucleus,
Propensities of receptors and genes activation and inactivation
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 rb(t) and rd(t), respectively, are equal for each copy:
The state of gene copy Gi (i = 1, 2) is Gi = 1 whenever NF-κB is bound to the promoter regulatory site, and Gi = 0 when the site is unoccupied. As a result the gene state G = G1 + G2 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 :
(1) At simulation time t; for given states of the A20 and IκBα genes, and number of active (bound) receptors B (M is the total number of receptors) we calculate the total propensity function r(t) of occurrence of any of the activation and inactivation reactions
(2) We select two random numbers p1 and p2 from the uniform distribution on [0, 1].
(3) Using the fourth order MATLAB solver we evaluate the system of 15 ODEs accounting for fast reactions, until time t + τ such that
(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.
In this step we determine which one of 8 potentially possible reactions occurs at time t+τ using the inequality
where ri(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.
In Figs. 8 and 9 we present evolution of selected variables during the persistent TNFα stimulation at high (10 ng/ml) and low (0.1 ng/ml) dose. We show both a single cell evolution and the average over 100 cells. The average is used to validate the model based on population data (Figs 3, 4, 7 and Figs S1, S2, S3, S4 in Additional file 3).
Figure 8. Evolution of main variables during persistent, 10 hours long stimulation with high dose TNFα = 10 ng/ml. Red line – typical single cell simulation, black line – average over 100 cells. Panel A: Number of active receptors, Panel B: IKKKa level, Panel C: IKKa level, Panel D: Total IκBα, Panel E: Nuclear NF-κB, Panel F: IκBα gene activity, Panel G: A20 gene activity, Panel H: IκBα mRNA, Panel I: A20 mRNA, Panel J: A20 protein.
Figure 9. Evolution of main variables during persistent, 10 hours long stimulation with low dose TNFα = 0.1 ng/ml, see Fig. 8 for description. Note, that at low dose stimulation first IKKa peak is smaller than for high dose but the subsequent peaks are well pronounced. As a result subsequent NF-κB oscillations are higher for low than high dose.
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.
Model fitting and parameters
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:
1) Start from a reasonable set of parameters, which produces a correct steady state in the absence of TNFα signal.
2) Proceed with the signal initiated by TNF downstream the autoregulatory loops.
3) 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.
A20 (TNFAIP3), TNF alpha inducible protein 3; IκBα, inhibitor of κB α subunit; IKK, IκB kinase; IKKK, IKK kinase; IRAK1 interleukin-1 receptor-associated kinase 1; MEKK3, mitogen-activated protein kinase kinase kinase 3; NF-κB, nuclear factor-κB; RIP, receptor interacting protein; TAK1, transforming growth factor-β-activated kinase 1; TNFα, tumor necrosis factor α; TNFR, TNF receptor; TRAF2, tumor necrosis factor receptor associated factor 2.
The author(s) declares that there are no competing interests.
TL, ARB and MK discussed the working hypothesis and wrote the paper. TL constructed and fitted the mathematical model. PP and MK discussed the mathematical model. PP wrote MATLAB programs. KP performed final numerical simulations and prepared the figures. ARB performed the IKK kinase activity experiment. All authors read and approved the final manuscript.
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.).
J Theor Biol 1991, 53:181-194. Publisher Full Text
PloS Biology 2006, 4:309. Publisher Full Text
Immunological Review 2006, 210:171-186. Publisher Full Text
Nelson G, Paraoan L, Spiller DG, Wilde GJC, Browne AM, Djali PK, Unitt JF, Sullivan E, Floettmann E, White MRH: Multi-parameter analysis of the kinetics of NF-κB signaling and transcription in single living cells.
Nelson DE, Ihekwaba AEC, Elliot M, Johnson JR, Gibney CA, Foreman BE, Nelson G, See V, Horton CA, Spiller DG, Edwards SW, McDowell HP, Unitt JF, Sullivan E, Grimley R, Benson N, Broomhead D, Kell DB, White MRH: Oscillations in NF-κB signaling control the dynamics of gene expression.
Science 2005, 308:52b. Publisher Full Text
Progress in Biophysics & Molecular Biology 2004, 86:5-43. Publisher Full Text
Hohmann H-P, Remy R, Piischl B, van Loon APGM: Tumor necrosis factors-α and β bind to the same two types of tumor necrosis factor receptors and maximally activate the transcription factor NF-κB at low receptor occupancy and within minutes after receptor binding.
Liu P, Jamaluddin M, Li K, Garofalo RP, Casola A, Brasier AR: Retinoic acid-inducible gene I mediates early antiviral response and toll-like receptor 3 expression in respiratory syncytial virus-infected airway epithelial cells.
Liu P, Choudhary S, Jamaluddin M, Li K, Garofalo R, Casola A, Brasier AR: Respiratory syncytial virus activates interferon regulatory factor-3 in airway epithelial cells by upregulating a RIG-I-toll like receptor-3 pathway.
Arnold R, Konig B, Galatti H, Werchau H, Konig W: Cytokine (IL-8, IL-6, TNF-ac) and soluble TNF receptor-I release from human peripheral blood mononuclear cells after respiratory syncytial virus infection.
PLoS Biology 2005, 3:1925-1938. Publisher Full Text
Current Biol 2003, 13:315-320. Publisher Full Text
Blonska M, Shambharkar PB, Kobayashi M, Zhang D, Sakurai H, Su B, Lin X: TAK1 is recruited to the tumor necrosis factor-α (TNF-α) receptor 1 complex in a receptor-interacting protein (RIP)-dependent manner and cooperates with MEKK3 leading to NF-κB activation.
Shim J-H, Xiao C, Paschal AE, Bailey ST, Rao P, Hayden MS, Lee K-Y, Bussey C, Steckel M, Tanaka N, Yamada G, Akira S, Matsumoto K, Ghosh S: TAK1, but not TAB1 or TAB2, plays an essential role in multiple signaling pathways in vivo.
Genes and Development 2006, 19:2668-2681. Publisher Full Text
Wertz IE, O'Rourke KM, Zhou H, Eby M, Aravind L, Seshagiri S, Wu P, Wiesmann C, Baker R, Boone DL, Ma A, Koonin EV, Dixit VM: De-ubiquitination and ubiquitin ligase domains of A20 downregulate NF-kB signaling.
J Chem Phys 2002, 117:6959-6969. Publisher Full Text