This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.056754v1 90/3/725    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Lipniacki, T. Articles by Kimmel, M. Search for Related Content PubMed PubMed Citation Articles by Lipniacki, T. Articles by Kimmel, M.

Stochastic Regulation in Early Immune Response

Tomasz Lipniacki ^*  , Pawel Paszek , Allan R. Brasier , Bruce A. Luxon  and Marek Kimmel  ^¶

^* Institute of Fundamental Technological Research, Warsaw, Poland; Department of Statistics, Rice University, Houston, Texas; Bioinformatics Program, and Department of Internal Medicine, University of Texas Medical Branch, Galveston, Texas; and ^¶ Institute of Automation, Silesian Technical University, Gliwice, Poland

Correspondence: Address reprint requests to Tomasz Lipniacki, Dept. of Statistics, Rice University, 6100 Main St., MS-138, Houston, TX 77005. E-mail: tomek{at}rice.edu; or E-mail: tlipnia{at}ippt.gov.pl.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MOTIVATION, MODEL FORMULATION,... RESULTS DISCUSSION APPENDIX I: JUSTIFICATION OF... APPENDIX II: MODEL EQUATIONS APPENDIX III: MODEL PARAMETERS APPENDIX IV: MODEL VALIDATION ACKNOWLEDGEMENTS REFERENCES

 
Living cells may be considered noisy or stochastic biochemical reactors. In eukaryotic cells, in which the number of protein or mRNA molecules is relatively large, the stochastic effects originate primarily in regulation of gene activity. Transcriptional activity of a gene can be initiated by transactivator molecules binding to the specific regulatory site(s) in the target gene. The stochasticity of activator binding and dissociation is amplified by transcription and translation, since target gene activation results in a burst of mRNAs molecules, and each copy of mRNA then serves as a template for numerous protein molecules. In this article, we reformulate our model of the NF-B regulatory module to analyze a single cell regulation. Ordinary differential equations, used for description of fast reaction channels of processes involving a large number of molecules, are combined with a stochastic switch to account for the activity of the genes involved. The stochasticity in gene transcription causes simulated cells to exhibit large variability. Moreover, none of them behaves like an average cell. Although the average mRNA and protein levels remain constant before tumor necrosis factor (TNF) stimulation, and stabilize after a prolonged TNF stimulation, in any single cell these levels oscillate stochastically in the absence of TNF and keep oscillating under the prolonged TNF stimulation. However, in a short period of 90 min, most cells are synchronized by the TNF signal, and exhibit similar kinetics. We hypothesize that this synchronization is crucial for proper activation of early genes controlling inflammation. Our theoretical predictions of single cell kinetics are supported by recent experimental studies of oscillations in NF-B signaling made on single cells.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MOTIVATION, MODEL FORMULATION,... RESULTS DISCUSSION APPENDIX I: JUSTIFICATION OF... APPENDIX II: MODEL EQUATIONS APPENDIX III: MODEL PARAMETERS APPENDIX IV: MODEL VALIDATION ACKNOWLEDGEMENTS REFERENCES

 
Typically, in mammalian cells, there are tens or hundreds of mRNA molecules of a given species and thousands of corresponding protein molecules, but typically only two homologous gene copies. This implies the discrete regulation of mRNA transcription, governed by the stochastic events of binding and dissociation of transcription factors to the regulatory site (1–6). Stochastic effects due to gene activation and inactivation are amplified by mRNA synthesis and protein translation. As a result, single gene activation may result in production of thousands of protein molecules, influencing the cell kinetics much more strongly than production and degradation of single mRNA or protein molecules. It was often assumed that the transcription rate of an inducible gene is a function of the nuclear concentration of a transcription factor. Such an approach is justified when the binding and dissociation events are very frequent, which may be assumed for prokaryotes, but is not justified for higher organisms. On the other hand, the large number of mRNA and protein molecules allows for a deterministic reaction-rate approximation in description of the time evolution of mRNA and protein levels, an approximation which is much more computationally efficient than the discrete model. Therefore, processes such as mRNA translation, formation, and degradation of protein complexes, and catalytic and spontaneous degradation, which involve a large number of molecules, may be modeled by ordinary differential equations (ODEs). As a contrast in prokaryotes, in which the mRNA molecules are very unstable (half-life time is of order of 1 min) and much less abundant, the stochasticity of formation, degradation, and translation of single mRNAs is of great importance (7–11). As a result, in prokaryotes there is a competition between stochastic effects caused by gene activation and mRNA processing.

Stochasticity in gene activation considerably influences cell kinetics and potentially leads to large cell-to-cell variability in mRNA and protein levels. Since each cell reacts to its own mRNA and protein levels, and not to the average levels in the population, the information about cell-to-cell variability could be crucial in the analysis of kinetics of single cell events like apoptosis. This is why, in this article on regulation of early immune response, we modify our previous simplified approach and model the transcriptional part of the regulatory network using a stochastic switch.

Nuclear factor B (NF-B) regulates numerous genes important for pathogen or cytokine inflammation, immune response, cell proliferation, and survival (reviewed in (12)). In mammals, the NF-B family of transcription factors contains five members but the ubiquitously expressed p50 and RelA heterodimer constitutes the most common inducible NF-B binding activity. In resting cells, p50-RelA heterodimers (referred to herein as NF-B) are sequestered in the cytoplasm by association with the members of another family of proteins called IB. This family includes several proteins but most of the IB-family inhibitory potential is carried by IB, whose synthesis is controlled by a highly NF-B-responsive promoter, generating autoregulation of NF-B signaling (13). Activation of NF-B requires degradation of IB, which allows NF-B to translocate into the nucleus, bind to B motifs present in promoters of numerous genes, and upregulate their transcription. NF-B activating signals converge on the cytoplasmic IB kinase (IKK), a multiprotein complex that phosphorylates IB, leading to its ubiquitination and then to its rapid degradation by the proteasome (reviewed in (14)). Activation of IKK kinase is induced by various extracellular signals including tumor necrosis factor- (TNF) and interleukin-1 (IL-1) through complicated, not fully resolved, transduction pathways. IKK inactivation is controlled by the zinc finger protein termed A20, which, like IB, is strongly NF-B responsive and generates a second autoregulatory loop in NF-B signaling (15). Mice deficient in A20 develop severe inflammation and cachexia, are hypersensitive to TNF, and die prematurely (16). 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 co-activators, pre-initiation factors, and RNA polymerase II, the mode of regulation for A20 and IB inhibitors appears to be distinct. Interestingly, chromatin immunoprecipitation assays have shown that the A20 and IB promoters are already bound by general transcription factors, co-activators, and RNA Pol II, waiting for the presence of NF-B binding to activate them (17,18). In these promoters, RNA polymerase II is stalled in an activated state; upon NF-B binding, RNA polymerase enters a competent elongation mode and thus is able to rapidly respond to NF-B binding. In this manner, these inhibitor members of the NF-B feedback loop are rapidly poised to respond to the presence of nuclear NF-B.

The first attempt to mathematically analyze the IB-NF-B signaling module was a one-feedback loop model that concentrated on the interplay between three IB isoforms (19). In contrast, our two feedback loop model (20) incorporates inhibitors operating at two levels in the signaling pathway: IB—the single species from IB family whose knockout is lethal (21), and which binds the majority of cytoplasmic RelA; and A20—which inhibits activity of IKK whose action leads to the second negative feedback loop in NF-B regulation. Incorporating the second inhibitor, A20, this latter model accurately predicts the profile of IKK activity.

In this article, we apply a stochastic switch together with deterministic reaction rate description by ODEs to model the NF-B regulatory network at the single cell level.

	MOTIVATION, MODEL FORMULATION, AND NUMERICAL IMPLEMENTATION

TOP ABSTRACT INTRODUCTION MOTIVATION, MODEL FORMULATION,... RESULTS DISCUSSION APPENDIX I: JUSTIFICATION OF... APPENDIX II: MODEL EQUATIONS APPENDIX III: MODEL PARAMETERS APPENDIX IV: MODEL VALIDATION ACKNOWLEDGEMENTS REFERENCES

 
In the last few years, several important studies of NF-B regulation at the single cell level have been performed (22–25). Using fluorescently tagged RelA and IB proteins, these experiments enable observation of intercompartment translocations of both proteins, showing a large heterogeneity in kinetics of cell responses to TNF or IL-1 stimulation. Nelson et al. (23) showed that the rate of nuclear RelA accumulation in response to TNF stimulation depends on the initial RelA/IB ratio. In cells with enhanced initial concentration of IB, the nuclear import is considerably slower. In accordance with Nelson et al. (23), Schooley et al. (25) found a broad distribution of nuclear RelA levels in the cell population at 10–90 min after IL-1 stimulation. Although these cited experiments indicate a large variability in cell kinetics, it is not straightforward to determine to which extent this variability is caused by stochastic regulation of gene expression, rather than being introduced by cell-to-cell variation in the amount of transfected plasmids (or their expression levels). In a recent experiment, Nelson et al. (26) proved the existence of long persisting single cell oscillations in NF-B signaling. Since cell-to-cell synchronization decreases with time, these oscillations either do not appear or appear strongly damped when observed at the population level, but as concluded by the authors they are important in control of expression of numerous NF-B-inducible genes. The point is that NF-B must circulate between nucleus and cytoplasm where it may become activated by a post-translational modification, such as phosphorylation. Inhibition of its nuclear export results both in NF-B dephosphorylation at Ser-536 after 3 h and in transitory B-dependent luciferase reporter gene expression that peaked after 5 h (26).

The model involves two-compartment kinetics of the activators IKK and NF-B, the inhibitors A20 and IB, and their complexes (Fig. 1). It is biologically justified to assume that the cytoplasmic complex IKK may exist in one of three forms (20):

1. Neutral (denoted by IKKn), synthesized de novo and specific to resting cells without an extracellular stimulus like TNF or IL-1.
   
2. Active (denoted by IKKa), arising from IKKn upon the TNF or IL-1 stimulation.
   
3. Inactive, but different from the neutral form, arising from IKKa possibly due to over-phosphorylation (denoted by IKKi).
   

View larger version (20K): [in this window] [in a new window]   FIGURE 1  Two-component feedback model of NF- B regulatory module. (A) Activation of NF-B module. NF-B is found in an inactivated complex in the cytoplasm with its inhibitor IB. TNF ligand acts as a switch to convert neutral IKK (IKKn) to the activated form, IKKa. IKKa, in turn, phosphorylates and degrades both free and NF-B complexed IB. Liberated NF-B enters the nucleus to induce transcription of IB and A20 genes. (B) NF-B-IB autoregulatory loop. The IB protein is rapidly resynthesized, enters the nucleus, and recaptures NF-B back into the cytoplasm. In the continued presence of IKKa, however, the resynthesized IB is continuously degraded, resulting in continued nuclear NF-B translocation. (C) The NF-B-A20 autoregulatory loop. A second level of negative autoregulation occurs with the resynthesis of A20. A20 is a ubiquitin ligase that degrades signaling intermediates coupling the TNF receptor with continuous IKK activation and directly associates with IKKa, converting it to catalytically inactive IKKi.

The transcriptional regulation of A20 and IB genes is governed by the same rapid elongation regulatory mechanism with rapid coupling between NF-B binding and transcription. The mechanisms for NF-B-dependent regulation of IB and A20 are based on control of transcriptional elongation, described earlier. 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. In fact, our experimental analysis of the kinetics of IB and A20 gene expression indicates that the patterns of mRNA expression are tightly coupled with NF-B presence in the nucleus, without appreciable time delay (28). This model cannot be extended to other NF-B-dependent genes that show distinct kinetics of induction. In this situation, the so-called late genes, like Naf-1 or NF-B2, show a peak of induction 6 h after NF-B binding. Activation of the late genes apparently requires the activation or binding of other rate-limiting regulatory factors (6). We assume that all cells are diploid, and both A20 and IB genes have two potentially active homologous copies, each of which is independently activated due to binding of the NF-B molecule to a specific regulatory site in the gene promoter. Following the literature (1–5,29), we made the simplifying assumption that each gene copy may exist only in two states, active and inactive. When the copy is active, the transcription is initiated at a high rate; and when the copy is inactive, transcription is inhibited, or it is initiated with some low basal rate. The gene copy inactivates when the NF-B molecule is removed from its regulatory site. We expect that this is due to the action of the IB molecules that bind to DNA-associated NF-B, exporting it out from the nucleus, rather than to thermal dissociation. The last would be very sensitive to temperature. The assumed mechanism of gene activation and inactivation is simplified, but it is also limited by our current knowledge of the subject. Fortunately, the details of activation mechanisms are not crucial for the model. The most important is the single assumption that transcription of A20 and IB genes turns on and off with probabilities determined by regulatory factors. This distinguishes the applied approach from the deterministic one in which the transcription speed is a function of concentrations of these factors.

The stochastic simulation algorithm (30) is an essentially exact procedure for simulation of the time evolution of processes like the ones described above. This exact algorithm becomes computationally intensive when the number of reacting molecules is large as in our case, i.e., 10⁵. In the last few years, several methods were proposed to accelerate the Gillespie algorithm. One is the -leap method (31,32), an approximate acceleration procedure, in which time is divided into -long intervals. It is required that is short enough so that the propensity functions for all reactions remain almost unchanged. Then, assuming that this condition is satisfied, the number of reactions in each reaction channel is a Poisson random variable, with parameter equal to the -interval propensity for the reaction channel considered. Another method is devised for the case when it is possible to separate system into slow and fast reaction channels (33,34). The idea is to approximate fast reactions either deterministically or as chemical Langevin equations, and to treat the slow reactions as stochastic events with time-varying reaction rates. Recently this idea has been improved by Cao et al. (35), who introduce a virtual fast system, being Markovian, which makes its analysis simpler. This improves the analysis in the cases when fast reactions are described by Langevin equations. In the case in which the fast reactions are fast enough to be described by the deterministic-rate equations, the method of Cao et al. (35) is reduced to that of Haseltine and Rawlings (33).

In this work, as in our recent article (29), we follow the method proposed by Haseltine and Rawlings (33) and split the reaction channels into fast and slow. Namely, we consider all reactions involving mRNA and protein molecules as fast and the reactions of gene activation and inactivation as slow. Fast reactions are approximated by the deterministic reaction-rate equations. The justification of this approximation will be given in the Appendices, where we compare it with the exact stochastic simulation algorithm in a simple example involving gene activation and inactivation, mRNA transcription, and protein translation. In this example, we will use the same reaction rates as appearing in our model.

According to the above, the mathematical representation of the model (see the Appendices) consists of 14 ordinary differential equations (ODEs) accounting for formation of complexes and their degradation, transport between nucleus and cytoplasm, transcription and translation, and of four algebraic equations accounting for propensity functions of binding and dissociation of NF-B molecules to regulatory sites in A20 and IB promoters. We assume that both A20 and IB genes have two homologous copies independently activated due to NF-B binding. It is assumed that in an infinitesimal time interval t, the probability P^b of NF-B binding to regulatory sites in each allele is proportional to the nuclear amount of NF-B,

(1)

NF-B dissociation probability, P^d, is a sum of a constant term and a term proportional to nuclear concentration of IB, which is capable of removing NF-B from regulatory sites in both A20 and IB genes,

(2)

It will be assumed that NF-B binding and dissociation proceed independently in each homologous gene copy and that binding and dissociation propensities r^b(t) = P^b(t, t)/t and r^d(t) = P^d(t, t)/t are equal for each copy. The state of gene copy Gⁱ (i = 1, 2) is Gⁱ = 1 whenever NF-B is bound to the promoter regulatory site, and Gⁱ = 0 when the site is unoccupied. As a result, the gene state G = G¹ + G² can be equal to 0, 1, or 2. In this approximation, the stochasticity of single cell kinetics solely results from discrete regulation of transcription of A20 and IB genes.

We assume that the transcription rate T_rate is the sum of a steady term and of an inducible term,

(3)

Note that Eq. 3 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 = 2 for most of the time. In such a case, the transcription would proceed at a maximum rate, c₂ + 2c₁. This rate has been measured for ß-actin by single RNA transcript visualization as 4 mRNA molecules per minute per one allele (2). In our calculations we assumed c₁ = 0.075 mRNA/s, which corresponds to 4.5 mRNA molecules per minute. It is worth noting that the protein production rate is proportional to the product of transcription and translation efficiencies. Therefore having solely information about the amount of protein one may not determine these two coefficients. When fitting the model, we found that even if c₂ = 0, transcription is regulated satisfactorily, and therefore to minimize the number of free parameters we assume c₂ = 0. For the same reason we set q₀ = 0 in Eq. 2.

In model computations, the amounts of all the substrates are specified in the numbers of molecules. Since we use the ODEs to describe most of the model kinetics, amounts of molecules are not integer numbers, but since these numbers are, in most cases, >>1, such a description is reasonable.

The implemented numerical scheme follows that of Haseltine and Rawlings (33), namely:

1. At simulation time t, for given states and of the A20 and IB genes, we calculate the total propensity function r(t) of occurrence of any of the binding or dissociation reactions,
   

   (4)

   
   

2. We select two random numbers p₁ and p₂ from the uniform distribution on (0, 1).
   
3. Using the fourth-order MATLAB solver we evaluate the system of 14 ODEs accounting for fast reactions, until time t + such that
   

   (5)

   
   

4. There are eight potentially possible reactions (NF-B may bind/dissociate to/from any of two alleles of A20 and IB genes), but at any given time, only four of them have nonzero propensities. In this step, we determine which one of eight potentially possible reactions occurs at time t + using the inequality
   

   (6)

   where r_i(t + ), i = 1, ..., 8 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.
   

Because of a large number of undetermined parameters, model fitting is not straightforward. We decided to carry out the fit manually, rather than to try to digitize the blot data and then to apply one of the fitting engines available. The first reason is that such quantification is by no means unique; the second is that, when fitting, we have to take into account diverse, and usually not precise, information coming from different researchers as well as our own intuitive understanding of the process. We first performed the fitting for the deterministic model (20), then for the stochastic model, leaving most of the parameters unchanged. Fitting of the deterministic model is considerably faster because it is not necessary to average the outcome over many runs. We applied the following fitting method:

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 along the autoregulatory loops.
   
3. Iterate item 2 until the fit to all the data is satisfactory.
   

The TNF signal first causes transformation of IKKn into IKKa. IKKa catalyses degradation of cytoplasmic (IB|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 IB genes. After being translated, A20 facilitates transformation of IKKa into IKKi, while IB 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 (IB|NF-B) and IB 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 fit is found, we observe that the set of parameters chosen to fit the data is by far 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 model onto the rest of the pathway is determined by their amounts multiplied by undetermined coupling coefficients. Hence, once we have a good fit, we may obtain another one using a smaller coupling coefficient and by proportionately enlarging the absolute level of the component. There also exist single parameters, the value of which can be changed over a very broad range without a significant influence on the time behavior of any substrate or complex for which the data are at our disposal. Values of all parameters are listed in the Appendices.

	RESULTS

TOP ABSTRACT INTRODUCTION MOTIVATION, MODEL FORMULATION,... RESULTS DISCUSSION APPENDIX I: JUSTIFICATION OF... APPENDIX II: MODEL EQUATIONS APPENDIX III: MODEL PARAMETERS APPENDIX IV: MODEL VALIDATION ACKNOWLEDGEMENTS REFERENCES

 
Stochasticity of the model implies that simulations performed with the same model parameters and the same initial conditions, are different. Since each simulation corresponds to the evolution in a single cell, therefore only the outcome averaged over many simulations corresponds to experimental data obtained for a population of cells in culture. As it will be seen in the Appendices, with parameters fitted, the proposed model is able to faithfully reproduce time-behavior of all variables for which the data is available: nuclear NF-B, cytoplasmic IB, A20, and IB mRNA transcripts, and total IKK and IKK kinase catalytic activity (IKKa) in both wild-type and A20-deficient cells.

In this section, we discuss the solutions of the model using the parameter set fitted to the experimental data from the experiments of Lee et al. (16) and Hoffmann et al. (19) on wild-type and A20-deficient mouse fibroblast cells.

Let us focus first on resting cell regulation. In the absence of TNF signal, IKK remains in its neutral (inactive) form, IKKn. This implies that it may not phosphorylate IB, which itself degrades at a relatively slow rate due to natural degradation. Until the amount of IB exceeds that of NF-B, which is assumed to be constant, equal to 60,000 molecules (see Appendices), all NF-B remains in cytoplasmic complexes with IB. In Fig. 2, we may observe that, because at t = 0 the number of cytoplasmic complexes of IB exceed 60,000 molecules, almost no NF-B is in the nucleus. When, due to degradation, the amount of IB falls below 60,000, some NF-B (typically a small fraction of cytoplasmic protein) enters the nucleus—where it may bind in a stochastic way to regulatory sites of IB and A20 promoters. Binding to IB results in a burst of IB transcription followed by an increase in IB protein level. In turn, free IB enters the nucleus, takes NF-B out of its regulatory site, and leads almost all NF-B back into the cytoplasm. Binding to A20 promoter results in a burst of A20 transcription, but A20 may not terminate the NF-B binding. As a result, each time NF-B enters the nucleus it causes a burst of IB, but not necessarily of A20. This implies that even if the (conditional) binding probabilities to IB and A20 promoters are equal, NF-B binds more frequently to IB promoters. If, however, NF-B binds to the A20 promoter before it binds to IB, it may remain there longer (since in this case binding does not lead to feedback reaction), and therefore may cause a larger burst of mRNA molecules. Although NF-B may bind to both homologous gene copies, in a resting cell (in which there are relatively few NF-B molecules in the nucleus), both copies are seldom simultaneously active. Let us note that since all cells in the population (in the absence of an external signal like TNF) behave asynchronously, the average IB and A20 transcript and protein levels in the cell population remain constant (plot not shown). The relatively low transcript level observed for IB and A20 in the population of unstimulated cells is usually explained by the existence of some basal transcription, independent of NF-B. However, our model shows that, even if the basal transcription rate is zero, the average IB and A20 mRNA levels may be positive and constant in time for unstimulated cells.

View larger version (24K): [in this window] [in a new window]   FIGURE 2  Numerical solution corresponding to a single (wild-type) resting cell. The amounts of substrates are given in numbers of molecules. D shows total IB, i.e., the sum of free and NF-B bound substrate in both the cytoplasm and nucleus.

View larger version (30K): [in this window] [in a new window]   FIGURE 3  Numerical solution corresponding to a single wild-type cell. The persistent TNF activation starts at 1 h. The amounts of substrates are given in numbers of molecules. B shows total IB, i.e., the sum of free and NF-B-bound substrate in both the cytoplasm and the nucleus.

View larger version (38K): [in this window] [in a new window]   FIGURE 4  Model simulations and measurements of Nelson et al. (26). Four single cell simulations (or experiments) are represented by differently colored lines. (A) Model simulations: nuclear NF-B after a persistent TNF activation starting at t = 0. The total amount of NF-B was elevated to 120,000 for this simulation. (B–D) Nuclear to cytoplasmic RelA-DsRed fluorescence normalized to highest peak intensity. Cells were treated with continual 10 ng/ml TNF. (B) SK-N-AS cells co-expressing RelA-DsRed and IB-EGFP. (C) SK-N-AS cells expressing RelA-DsRed and EGFP control. (D) HeLa cells expressing both RelA-DsRed and IB-EGFP.

View larger version (29K): [in this window] [in a new window]   FIGURE 7  Numerical solution averaged over a population of 500 wild-type cells; all parameters as in Fig. 3. The amounts of substrates are given in numbers of molecules per cell.

View larger version (49K): [in this window] [in a new window]   FIGURE 5  Oscillations of nuclear NF-B after a persistent TNF treatment for six different levels of total NF-B: 10,000 / 20,000 / 30,000 / 60,000 / 120,000 / 180,000.

View larger version (35K): [in this window] [in a new window]   FIGURE 6  Correlation of the total NF-B level with the peak-to-peak timing of NF-B oscillations. (Left column) Nelson et al. (38) experiments on SK-N-As cells. (Right column) Model predictions. The line represents the average calculated based on 100 simulations performed at each point, corresponding to the total amount of NF-B, equal to six different levels: 10,000 / 20,000 / 30,000 / 60,000 / 120,000 / 180,000. The dots corresponds to 500 single simulations, each done for the amount of total NF-B randomly selected from uniform distribution on 10,000 and 180,000.

View larger version (20K): [in this window] [in a new window]   FIGURE 9  Scatter plots of total A20 versus nuclear NF-B at six times after start of TNF stimulation.

In a cell co-expressing IB-EGFP and RelA-DsRed (Fig. 4 B), one can observe large differences in duration of the first peak. As stated by Nelson et al. (26), these differences are due to the sensitivity of the system to variations in IB-expression efficiency. Cells with higher expression of IB-EGFP show statistically prolonged NF-B nuclear activity followed by a longer period in which NF-B is absent in the nucleus. We do not observe, however, such strong sensitivity of the first peak duration to the IB expression efficiency in our model (data not shown). It is possible that the strong extension of the first peak duration and generally more-complex pattern of peak amplitude is because the kinetics of IB-EGFP is slower than the kinetics of endogenous IB. One might expect that the fusion protein is degraded at a lower rate (under the TNF stimulation) and it is not as strongly induced by NF-B. In Fig. 1 (26), we can observe that in cells not expressing IB-EGFP (initially pure red), the endogenous IB is degraded in the first 6 min (and as a result RelA-DsRed enters the nucleus), but in cells expressing IB-EGFP, it remains fairly abundant even after 60 min. Similarly, in cells ((26); see Fig. 6S, in that article's Supplementary Data) which, at t = 0, co-express IB-EGFP and RelA-DsRed, the endogenous IB rebuilds more quickly than IB-EGFP. At 300 min, RelA-DsRed is taken out from the nucleus—suggesting that the endogenous IB is rebuilt, but that the IB-EGFP is still absent (cells are purely red).

In Fig. 7, we show the same kinetics as in Fig. 3, but averaged over 500 simulation runs or over a population of 500 identical cells. The averaged kinetics is considerably different from kinetics of a particular cell. Despite the fact that each cell oscillates under extended TNF stimulation, all quantities averaged over many cells converge to steady-state levels. This effect is caused by a growing desynchronization of cells. Before the TNF signal, the cells are desynchronized (the single cell simulations are started at 20–30 h before the TNF signal), and it is the TNF signal that induces cell synchronization. The sharp decrease in total IB, followed by buildup of nuclear NF-B, is observed in all simulations (Fig. 8). Thereafter, the stochastic nature of NF-B binding causes that the peaks of gene activity to not-match, and this desynchronizes cell kinetics. The activity of the IB and A20 genes, averaged over a relatively large population, corresponds much better to the nuclear NF-B level than it does in the single cell case.

View larger version (17K): [in this window] [in a new window]   FIGURE 8  Scatter plots of total IB versus nuclear NF-B at six times after start of TNF stimulation.

The scatter plots in Fig. 10 show abundance of I B and A20 proteins compared to their mRNA at three time points. We observe that initially there are relatively few (<100) IB mRNA molecules per cell, although the IB protein (which is mostly complexed with NF-B) is abundant. Then at 30 min most of the IB protein is degraded, but the number of message molecules is large due to NF-B-induced transcription. The kinetics of A20 transcript and protein is different. Initially there is a small amount of both protein and mRNA, then a growing amount of transcript is followed by the growing amount of protein. We observe broadening of the distribution in time, caused by desynchronization of cells due to stochasticity. Despite the fact that the scatter plots presented in Figs. 7–9 reveal that variability among cells grows in time, one may expect that there is a sizable fraction of cells, the evolution of which is close to the averaged evolution. In fact, however, none of the cells behave like the average. In Fig. 11, we show that single cell trajectories keep oscillating when the equilibrium distribution is reached, whereas the average trajectory resulting from averaging over 500 simulated cells stabilizes.

View larger version (22K): [in this window] [in a new window]   FIGURE 10  (Scatter plots, left) Total IB protein versus IB message. (Scatter plots, right) A20 protein versus A20 message at three times after start of TNF stimulation.

View larger version (34K): [in this window] [in a new window]   FIGURE 11  Trajectories projected on the (total IB, nuclear NF-B, time) hyperplane. (Thin lines) Three single cell trajectories; (thick line) average trajectory.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MOTIVATION, MODEL FORMULATION,... RESULTS DISCUSSION APPENDIX I: JUSTIFICATION OF... APPENDIX II: MODEL EQUATIONS APPENDIX III: MODEL PARAMETERS APPENDIX IV: MODEL VALIDATION ACKNOWLEDGEMENTS REFERENCES

The method of modeling molecular pathways applied by us follows that of Haseltine and Rawlings (33) and splits the reaction channels into fast and slow. It combines deterministic reaction-rate description (ODEs) for fast reaction channels (in this case, processes involving large number of molecules) with a stochastic switch for gene activity. This approach is computationally much more efficient than the stochastic simulation algorithm (30), but still accurate enough for our purposes (see Fig. 12).

View larger version (24K): [in this window] [in a new window]   FIGURE 12  Exact stochastic simulation (thin line) compared with stochastic-deterministic approximation (thick line). (A) mRNA profile; (B) protein profile. The transcription, translation, and degradation rates are the same as for A20 gene; see Table 2.

View this table: [in this window] [in a new window]   TABLE 2  A20 and IB synthesis and degradation, IKK dynamics, and total amount of free and complexed NF-B