# Statistical ensemble analysis for simulating extrinsic noise-driven response in NF-κB signaling networks

Jaewok Joo, Steven J. Plimpton, Jean-Loup Faulon

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jaewok Joo, Steven J. Plimpton, Jean-Loup Faulon. Statistical ensemble analysis for simulating extrinsic noise-driven response in NF-κB signaling networks. BMC Systems Biology, 2013, 7, 10.1186/1752-0509-7-45. doi: 10.1186/1752-0509-7-45

HAL Id: hal-01601848

https://hal.science/hal-01601848

Submitted on 29 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Statistical ensemble analysis for simulating extrinsic noise-driven response in NF-κB signaling networks

Jaewook Joo¹*, Steven J Plimpton² and Jean-Loup Faulon³

**Abstract**

**Background:** Gene expression profiles and protein dynamics in single cells have a large cell-to-cell variability due to intracellular noise. Intracellular fluctuations originate from two sources: *intrinsic* noise due to the probabilistic nature of biochemical reactions and *extrinsic* noise due to randomized interactions of the cell with other cellular systems or its environment. Presently, there is no systematic parameterization and modeling scheme to simulate cellular response at the single cell level in the presence of extrinsic noise.

**Results:** In this paper, we propose a novel statistical ensemble method to simulate the distribution of heterogeneous cellular responses in single cells. We capture the effects of extrinsic noise by randomizing values of the model parameters. In this context, a statistical ensemble is a large number of system replicates, each with randomly sampled model parameters from biologically feasible intervals. We apply this statistical ensemble approach to the well-studied NF-κB signaling system. We predict several characteristic dynamic features of NF-κB response distributions; one of them is the dosage-dependent distribution of the first translocation time of NF-κB.

**Conclusion:** The distributions of heterogeneous cellular responses that our statistical ensemble formulation generates reveal the effect of different cellular conditions, e.g., effects due to wild type versus mutant cells or between different dosages of external stimulants. Distributions generated in the presence of extrinsic noise yield valuable insight into underlying regulatory mechanisms, which are sometimes otherwise hidden.

**Keywords:** Statistical ensemble, Extrinsic noise, Cell to cell variability, NF-κB signal transduction network

**Background**

Single cell imaging generated a surge of interest in the intracellular dynamics of biochemical species, uncovering significant cell-to-cell variations in gene expression [1-8] and protein dynamics [9,10]. This variability originates from intrinsic [1-8] and extrinsic noise [3,6,10] and critically affects cellular decision-making processes [9-13]. Moreover, cellular response averaged over a population of cells is oftentimes noticeably different from the responses of single cells. The variability in the latter contains rich information regarding the regulatory mechanisms in operation. Here, we present a novel computational method to predict the distribution of extrinsic noise-driven heterogeneous cellular

responses and to unravel discrepancies between single-cell versus population-averaged responses.

Both intrinsic and extrinsic noise are the source of the large cell-to-cell variability in cellular responses [14]. *Intrinsic* noise refers to the pure probabilistic nature of individual biochemical reactions occurring within a cell. When the number of intracellular constituents is large, the cell's behavior is well approximated by its expectation value according to the law of large numbers. But at the single-cell level, the number of molecules of certain species critical to a particular biochemical pathway can be small, and the range of statistical variation in the system needs to be considered [1-8]. *Extrinsic* noise refers to random interactions of the cell with other cells or its environment. Extrinsic fluctuations can originate from cells undergoing different stages of their cell cycle [15], fluctuations in the number of transcriptional regulators upstream of the signaling pathway of interest [3,6,9,10],

<footnote>* Correspondence: jjoo1@utk.edu
¹Department of Physics and Astronomy, University of Tennessee, Knoxville 37996, USA</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
and cell-to-cell variability in the copy number of proteins inherited from parent cells during cell division [10]. Extrinsic noise can affect the dynamics of cellular constituents locally in a specific signaling pathway or globally over the entire cell. In Figure 1, we summarize the effects of intrinsic and extrinsic fluctuations in the NF-κB signaling networks. The full effect of extrinsic noise should include “all” external stochastic effects that influence the cell, particularly the temporal fluctuations in the cellular kinetic conditions. However, in Ref. [10], Spencer et al. identified the most important source of extrinsic noise as the protein copy number inherited from the parent cell during cell division. Large cell-to-cell variations in the copy number of enzyme and regulatory protein could randomize the likelihood and the speed of any intracellular biochemical reaction. This means we can effectively “lump” all the effects of protein copy number variations into variations in kinetic rate constants. This is an attractive approach, because rate constants are an input into a variety of biochemical pathway modeling techniques.

A pathway modeling framework that uses deterministic or stochastic differential equation models requires a priori knowledge of the structure of the biochemical reaction network, mathematical functional forms for the biochemical reactions, and associated reaction rate constants. Since limited or incomplete information is often all that is available to modelers, a computational model is often parameterized by using a nonlinear fitting algorithm. A conventional parameterization scheme identifies a single set of kinetic parameter values by minimizing the χ² distance between experimental data and a prediction made by the model. Sloppy Cell and other similar parameterization algorithms include experimental errors in the parameterization by fitting to a rather large experimental error bar [16]. But both

conventional and Sloppy Cell parameterization schemes assume a deterministic and homogeneous biological response to a stimulus and aren't designed to handle the heterogeneous, stochastic behavior of single cells and its dependence on extrinsic noise.

In order to capture extrinsic noise and its effect on intracellular response, we propose a novel parameterization method, the “statistical ensemble” (SE) scheme, named after a key concept in statistical physics [17]. A cell is regarded as a complex system comprising a large number of components and elementary interactions among them. A population of cells consists of a large number of replicates, each with different microscopic intracellular states. The statistical ensemble average, or macroscopic observable, is equated with the cellular response averaged over the population of cells. The ensemble is generated by assigning randomly sampled values of kinetic rate constants and copy numbers of regulatory proteins to each cell in the ensemble. All other external noisy systems that interact with the cell, but which are not modeled explicitly, are treated as extrinsic noise. The effect of the noise is included in the sampling that produces the randomized microscopic state of each cell in the ensemble.

A key point is that the resulting dynamic response of the ensemble of cells is no longer a single output but is a distribution of heterogeneous responses. Each response can be computed independently, which allows for parallelism in the computation. An equal weight is assigned to the response from each replicate, to calculate the ensemble averaged cellular response. The SE scheme thus enables modeling of the irregular, dissimilar, and diverse individual cellular behaviors while reproducing the macroscopically observable population-level response.

In the most general sense, the success of the SE scheme depends on identifying and characterizing the biologically

<figure><img src="image_3.png" /><figcaption>**Figure 1** Intrinsic and extrinsic noise as the source of the cell-to-cell variability in cellular responses in the NF-κB signaling networks.</figcaption></figure>

Intrinsic noise refers to the pure probabilistic nature of individual biochemical reactions in the signaling networks. Extrinsic noise refers to random interactions of the signaling networks with the external stochastic systems and originates from three sources: (i) fluctuating number of transcriptional regulators upstream of the signaling networks, (ii) fluctuating number of proteins inherited from parent cells, and (iii) different stages of their cell cycle.
correct distribution of extrinsic noise for the system of interest, so that its effects can be encoded in the random sampling of cellular microstates. In this work, we use experimental population-level data to parameterize the range of feasible kinetic rate constants and copy numbers of specific molecules, and then sample uniformly around the mid-point of the range to generate cellular microstates. We illustrate the power of even this simplified SE approach for modeling the NF-κB signaling system.

NF-κB is a pleiotropic regulator of gene control and plays significant roles in various cellular functions such as differentiation of immune cells, development of lymphoid organs, and immune activation [18-20]. NF-κB shuttling between the nucleus and cytoplasm is auto-regulated by the NF-κB signaling module, which consists of IκB (inhibitor κB), IKK (IκB kinase), and NF-κB. In the absence of stimulus, IκB forms a hetero-dimeric complex with NF-κB, preventing NF-κB from entering into the nucleus. Upon stimulation, phosphorylated IKK catalyses the degradation of IκB from the IκB-NF-κB complex and frees up NF-κB whose nuclear localization initiates transcription of NF-κB target genes such as inflammatory cytokines (TNFα, IL-1, IL-6), chemotactic cytokines (MIP-1α), Th1 and Th2 response activation (IFN and IL-10), and lastly, but most importantly, negative regulators (IκBα, IκBβ, IκBε, and A20) which terminate the NF-κB signaling. Based on current knowledge of NF-κB signaling, Hoffmann et al. constructed a complex biochemical reaction network for the NF-κB signaling pathway consisting of IKK, NF-κB, and three IκB isoforms and transformed it into a set of ordinary differential equations with dozens of unknown kinetic parameter values [21]. After identifying a single set of parameter values yielding the best fit of the model to population level experimental data, they used their model to elucidate the role of each of three IκB isoforms: IκBα induces oscillatory shuttling of NF-κB while IκBβ and IκBε damp the oscillations [21]. Lipniacki et al. extended the model, showing that an additional negative regulator A20 has a definitive role as a NF-κB signal terminator, by deactivating IKK phosphorylation [22-24]. Using fluorescence microscopy, Nelson et al. and several other groups showed a remarkably heterogeneous intracellular response for this signaling network at the single-cell level; some cells exhibited sustained oscillatory shuttling of NF-κB while others exhibited non-oscillatory behavior [25-33].

In this paper, we model extrinsic noise via randomization of the kinetic parameters of the IKK-NF-κB-IκB-A20 signaling system and predict several distributions of dynamical NF-κB responses. The signaling network we model is shown in Figure 2 and consists of IKK, cytoplasmic and nuclear NF-κB, and two groups of negative regulators (three isoforms of IκB and A20). Using the statistical ensemble (SE) scheme, we demonstrate that extrinsic noise, modeled as fluctuations in kinetic parameter values, can generate the

observed experimental population-level response as the SE average, as well as a heterogeneous distribution of individual cellular responses. In section Results.A we show that the SE average of key biochemical species concentrations in the NF-κB signaling network can be accurately fit to experimental population-level data for wild type and various mutant cases. In section Results.B, we predict the distributions of various dynamic characteristics of NF-κB cellular responses. In section Results.C, we make a prediction about dosage-dependent NF-κB responses in single cells, i.e., the dosage-dependent distribution of various NF-κB dynamic characteristics in individual cells. Lastly, in section Results.D, we predict that both dose-response curves from individual cells and their SE average are sigmoidally shaped.

## Results

**Statistical ensemble average of key biochemical species concentrations in the NF-κB signaling network is fit to experimental population-level data**

### The wild type case

For this reaction pathway, the statistical ensemble (SE) scheme generates significant cell-to-cell variability in protein dynamics. Yet the SE averages agree well with population-level experimental data (Electro Mobility Shift Assay (EMSA) or western blot) for key biochemical species concentrations as shown in Figure 3. For the nuclear NF-κB profiles in Figure 3(A), the first translocation times (timing of the first peak) of the individual NF-κB profiles (in blue) are almost identical, while the first maxima (amplitude of the first peak) vary significantly with

<figure><img src="image_4.png" /><figcaption>**Figure 2 Biochemical network model for the IKK-IκB-NF-κB-A20 signaling module.** Top panel: A schematic description of our comprehensive model for NF-κB signaling. The arrows indicate activation and the perpendicular lines denote inhibition. Bottom panel: the model consists of IKK (IκB kinase), IκB isoforms (IκBᵢ, i = α, β, ε), and A20. NF-κBn and IκB/n denote their nuclear components. Squares are for proteins; hexagons are for mRNA. Black arrows indicate either association or dissociation or degradation of proteins; red (blue) arrows denote mRNA (protein) synthesis.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3** Individual time-series curves (blue lines) and the ensemble average (red line) of key protein concentration for an ensemble of 1000 replicates of the wild type NF-κB signaling system. Computational results are compared side-by-side with population-level experimental data from Ref. [20]. Panel (A): nuclear concentration of NF-κB. Panels (B), (C), (D) are respectively cytoplasmic concentrations of IκBα, IκBβ, and IκBε proteins.</figcaption></figure>

a variance up to 100% of the SE average (in red). However, both the timings and amplitudes of subsequent peaks exhibit significant cell-to-cell variability. Consequently the SE average is a strongly damped oscillatory pattern with rapid decay of subsequent peak amplitudes. Thus, the effect of extrinsic noise on this observable is a “masking effect of averaging over a population of asynchronous curves”, just as for intrinsic noise [34]. The large variation in the first-peak amplitude of nuclear NF-κB concentration in Figure 3(A) originates from the IKK profile in Figure 4(C), where the IKK concentration time courses from individual cells also exhibit significant differences in their first maximum. This induces large variation in the first minimum of IκB isoforms as shown in Figures 3(B)-(D). Thus, the cell-to-cell variation in kinetic rate constants regulating the levels of both pre-activated IKK (IKKn) and activated IKK

(IKKa) is the source of similar variation in the first maximum of nuclear NF-κB concentration [35]. Likewise, the asynchronous behavior of the individual nuclear NF-κB profiles after two hours, as shown in Figure 3(A), originates from the cell-to-cell variability in the second-peak amplitude of the IκB isoforms in Figures 3(B)-(D).

### The mutant case - double knocked-out IκB isoforms and knocked-out A20

To simulate the dynamics of mutants, we set the mRNA synthesis rates for two of the three IκB isoforms and A20 to zero. For the IκBβ and IκBε knocked-out mutant shown in Figure 5(A), the peaks of the SE average correspond closely to the peaks of population-level experimental data (EMSA) at 15 min, 2.5 hours, 4 hours, and 5.5 hours. The individual profiles of nuclear NF-κB concentration are
<figure><img src="image_6.png" /><figcaption>**Figure 4** Individual time-series curves (blue lines) and the ensemble average (red line) of key protein concentrations are obtained for an ensemble of 1000 replicates of a wild type (A, C) and an A20 knocked-out mutant (B, D). Computational results are compared with population-level experimental data from Ref. [36]. Top panels: nuclear concentration of NF-κB. Bottom panel: IKK concentration.</figcaption></figure>

much more oscillatory (about half the curves exhibit sustained oscillations as shown in Figure 6) than for the wild type data (only 10% are sustained oscillations in Figure 6). But, the SE average of this mutant is a damped oscillatory pattern, with a bit more dynamic variation than that of the wild type. This is again mainly due to “the masking effect of averaging over a population of asynchronous curves”. For the IκBα and IκBε knocked-out and the IκBα and IκBβ knocked-out mutants shown in Figure 5(B) and 5(C), the SE averages of nuclear NF-κB show a “single-peaked” pattern similar to the population-level EMSA data, though the timings of the peaks differ by 1 hour. The single-peak amplitudes vary significantly with a variance as large as 100% of the SE average. For the A20 knocked-out mutant in Figure 4(B) and 4(D), both the SE averages of nuclear NF-κB and IKK profiles exhibit single-peaked patterns in good agreement with the population-level experimental data. Again, the individual nuclear NF-κB profiles differ significantly. For all the mutants, though their SE averages for nuclear NF-κB

profiles exhibit simple dynamic patterns, the cell-to-cell variability is large due when extrinsic noise is included in the model.

### Dependence of SE average on heterogeneity

In Figure 7 we show how to use population-level experimental data as a constraint when choosing a heterogeneity factor $\chi$, defined as the interval size of the uniform distribution from which kinetic rate constants are sampled, as inputs to the pathway model. Centering the kinetic rate constants at their reference values, we vary $\chi$ and observe how heterogeneous the individual cell profiles of nuclear NF-κB become. Note that the SE average of nuclear NF-κB becomes less oscillatory for higher values of $\chi$ in Figure 7. For a small $\chi$ = 10% in Figure 7(A), all individual curves remain in phase with each other, making the SE average also highly oscillatory. For higher values of $\chi$ = 50% and $\chi$ = 70% in Figures 7(C) ($\chi$ = 50%) and 7(d) ($\chi$ = 70%), a large fraction of individual curves are sustained oscillations, but quickly
<figure><img src="image_7.png" /><figcaption>**Figure 5** Individual time-series curves (blue lines) and the ensemble average (red line) of key protein concentrations for an ensemble of 1000 replicates of a *IκB* double gene *knocked-out* mutant. Computational results (left column) are compared with population-level experimental data (right column) from Ref. [20]. Panel (A): *IκBβ* and *IκBε* knocked-out mutant. Panel (B): *IκBα* and *IκBβ* knocked-out mutant. Panel (C): *IκBα* and *IκBε* knocked-out mutant.</figcaption></figure>

become out of phase, resulting in an SE average that is strongly under-damped. Because higher χ values cover a larger sampling space, individual nuclear NF-κB curves bifurcate into different classes of patterns: some are sustained oscillatory while others are single-peaked. Thus, if the population-level experimental data exhibit sustained oscillations versus damped oscillations versus single-peak profiles, the variation in single-cell profiles induced by χ can be used to guide sampling from an appropriate range of heterogeneity when generating input rate constants.

In this subsection, we showed how the SE method with its many replicates is a model for a population of cells in a heterogeneous set of intracellular states. By varying the

heterogeneity factor χ for sampling kinetic parameters used as inputs to the pathway model, fits to experimental data can be produced even when population-averaged data and single-cell data exhibit different characteristics, as in the NF-κB signaling system. In the next subsection, we discuss the distribution of single-cell NF-κB responses in more detail.

### Prediction of distributions of individual cellular responses for the wild type and mutants

#### Distributions of dynamic features

In Figure 8, we summarize the output of our SE computational model for distributions of single-cell responses for the wild type and the mutants discussed in the previous sub-section. Six dynamic features are shown: the
<figure><img src="image_8.png" /><figcaption>**Figure 6** Distributions of four dynamic patterns of the individual time-series curves of nuclear NF-κB profiles for an ensemble of 1,000 replicates of the wild type, A20 knocked-out mutant, and three IκB genes double knocked-out mutants. A few examples of four dynamic patterns are plotted in the top panel: (A) single-peaked pattern (blue), (B) under-damped oscillation (red), (C) hyperbolic pattern (black), and (D) sustained oscillation (yellow) where color within a parenthesis denotes color in the bottom panel. Individual time-series curves are classified as one of the four dynamic patterns.</figcaption></figure>

amplitude of the first peak (First Maximum), the timing of the first peak (First Translocation Time), the time between the first and the second peaks (First Period), the level of the first minimum (First Minimum), the amplitude of the second peak (Second Maximum), and the asymptotic Steady State value. Surprisingly, for each dynamic feature, there is a significant amount of overlap between the distributions for the wild type and those of the mutants. This implies that if we used a conventional modeling scheme which fits a single set of parameter values and outputs a single representative time-series of intracellular response, we could draw incorrect conclusions as to the effect of a

knocked-out gene on cellular response. To avoid this, we compute the entire distribution of responses and look for significant changes when genes are knocked out. In Figure 8(A), the distributions of the First Maximum are the same for both the mutants and the wild type. This dynamic feature is thus not an indicator of the physiological defects caused by the knock-out genes. In Figure 8(B), the distribution of the First Translocation is shifted to the right for the A20 knocked-out mutant and to the left for IκBβ and IκBε double knocked-out mutant, whereas the wild type and two other mutants have similar distribution. In Figure 8(C), only the wild type and the IκBβ
<figure><img src="image_9.png" /><figcaption>**Figure 7** Dependence of the individual time-series curves (blue lines) and the statistical ensemble average (red line) of nuclear NF-κB profiles for a mutant with IκBβ and IκBε genes double knocked-out, on the heterogeneity factor χ (the interval size of the uniform distribution or kinetic rate parameters). (A) χ = 10%; (B) χ = 30%; (C) χ = 50%; (D) χ = 70%.</figcaption></figure>

and IκBε double knocked-out mutant have well-defined periods of roughly 2 hours; the First period of other mutants is too broadly distributed to define an average. In Figure 8(D), the ratio of the First Minimum to the First Maximum indirectly measures the spikiness of the oscillations; the smaller the ratio, the spikier the temporal profile becomes. Only the wild type and the IκBβ and IκBε double knocked-out mutant exhibit a spiky response. In Figure 6(F), the ratio of the Steady State to the First Maximum provides useful information about the relative magnitude and strength of the negative regulators of IκB isoforms and A20. Since the distributions of the First Maximum are the same for the wild type and mutants, we conclude that the smaller steady-state level of nuclear NF-κB concentration infers stronger negative feedback. The mutants ordered by steady-state level are as follows: A20 knocked-out mutant < IκBα and IκBε knocked-out mutant < IκBα and IκBβ knocked-out mutant < IκBβ and IκBε knocked-out mutant < wild type. The relative strength of the negative regulators can then be inferred: A20 > IκBα > IκBε > IκBβ. Of course, this ordering is consistent with the choice of nominal values for the respective kinetic rate constants, as listed in Table 1.

### Distribution of dynamic patterns

The individual time-series of the nuclear NF-κB concentrations can be classified into one of four dynamic patterns

(damped oscillation, sustained oscillation, single peaked, and monotonic-increasing patterns) as shown in Figure 6. The underlying mechanism for each dynamic pattern is rather simple. The monotonic-increasing (or over-damped) pattern originates from strong negative feedback loops, while the single-peaked pattern results from weak negative feedback loops. The oscillatory patterns arise from intermediate-strength negative feedback loops. But it remains an open question to correlate each dynamic pattern with a specific cellular physiology [37-39]. To elucidate this connection, we stimulate the ensemble of NF-κB signaling networks with the same signal strength (TR = 1), for both the wild type and mutants. We then classify a thousand individual temporal profiles into one of the four dynamic patterns. The distributions of the patterns are represented by bar graphs in Figure 6 which shows that both the wild type and mutants exhibit at least two different patterns of response under the same strength of stimulation. For the wild type, most of the nuclear NF-κB profiles have a damped-oscillatory pattern, with less than 10% of the profiles as sustained-oscillatory. This indicates a damped-oscillatory response is the most probable, and it is robust against perturbation of the network parameter values. For the mutant with a knocked-out A20 gene, both single-peaked and damped-oscillatory patterns are nearly equally probable. But the damped oscillatory profiles are very similar to a single-peaked pattern. Thus for this mutant, a damped-
<figure><img src="image_10.png" /><figcaption>**Figure 8** Distributions of six dynamic features of nuclear NF-κB profiles, for an ensemble of 1,000 replicates of the wild type (black lines), A20 knocked-out mutant (red lines), *ΙκΒα* and *ΙκΒβ* genes double knocked-out mutant (blue lines), *ΙκΒα* and *ΙκΒε* genes double knocked-out mutant (yellow lines), and *ΙκΒβ* and *ΙκΒε* genes double knocked-out mutants (green lines). The six dynamic features are First Maximum (the amplitude of the first peak) in panel (A), First Translocation Time (the timing of the first peak) in panel (B), First Period (the time between the first and the second peaks) in panel (C), Ratio of First Minimum to First Maximum (ratio of the first minimum value to the first maximum value) in panel (D), Ratio of Second Maximum to First Maximum (ratio of the second peak amplitude to the first maximum value) in panel (E), and Ratio of Steady State to First Maximum (ratio of the steady state level to the first maximum value) in panel (F).</figcaption></figure>

oscillatory response occurs in a region of the parameter space where the negative regulators are not strong enough to induce the oscillatory pattern. For the mutant with *ΙκΒβ* and *ΙκΒε* genes double knocked-out, sustained-oscillatory and damped-oscillatory patterns are equally probable responses. The damped-oscillatory responses in this mutant are very different from those in the mutant with a knocked-out A20 gene, and are more similar to a sustained oscillation. The fraction of sustained-oscillatory responses (about 50%) dramatically increases in comparison to the wild type case (less than 10%). For mutants with *ΙκΒα* and *ΙκΒβ* genes double knocked-out and with *ΙκΒα* and *ΙκΒε* genes double knocked-out, their respective distributions are similar to that of the mutant with the A20 gene knocked-out. As shown in Figures 5(B) and 5(C) and Figure 4(B), both the

individual profiles and the statistical ensemble average of the nuclear NF-κB concentrations for all these mutants (A20 gene knocked-out, *ΙκΒα* and *ΙκΒβ* genes double knocked-out, *ΙκΒα* and *ΙκΒε* genes double knocked-out) exhibit similar single-peaked patterns. In summary, there are two distinctive groups exhibiting two respective patterns of nuclear NF-κB profile response: the first group, consisting of the wild type and the *ΙκΒβ* and *ΙκΒε* double knocked-out mutant, is dominated by highly oscillatory responses. The second group, consisting of the A20 knocked-out, the *ΙκΒα* and *ΙκΒβ* double knocked-out, and the *ΙκΒα* and *ΙκΒε* double knocked-out mutants, shows mostly single-peaked (non-oscillatory) responses.

In this subsection, we used the SE method to generate distributions of various dynamical responses from a large
<table><caption>Table 1 Biochemical reactions and their associated reaction rate constants in the computational model of the NF-κB signaling network</caption><thead><tr><th>Reactions</th><th>I</th><th>II</th><th>III</th><th>IV</th><th>V</th></tr></thead><tbody><tr><td>IKKa + IkBα → IKKa:IKBα</td><td>Aα</td><td>[a]</td><td>0.2</td><td>[1]</td><td>0.1813</td></tr><tr><td>IKKa + IkBβ → IKKa:IKBβ</td><td>Aβ</td><td>[a]</td><td>0.05</td><td>[3]</td><td>0.02997</td></tr><tr><td>IKKa + IkBε → IKKa:IKBε</td><td>Aε</td><td>[a]</td><td>0.05</td><td>[3]</td><td>0.04244</td></tr><tr><td>IKKa + IkBα:NF-κB → IKKa:IKBα:NF-κB</td><td>Bα</td><td>[a]</td><td>1</td><td>[1]</td><td>1.024</td></tr><tr><td>IKKa + IkBβ:NF-κB → IKKa:IKBβ:NF-κB</td><td>Bβ</td><td>[a]</td><td>0.25</td><td>[3]</td><td>0.3683</td></tr><tr><td>IKKa + IkBε:NF-κB → IKKa:IKBε:NF-κB</td><td>Bε</td><td>[a]</td><td>0.25</td><td>[3]</td><td>0.42</td></tr><tr><td>NF-κBn → NF-κBn + A20t</td><td>C1</td><td>[b]</td><td>0.0000005</td><td>[1]</td><td>0.000000506</td></tr><tr><td>0 → A20t</td><td>C2</td><td>[c]</td><td>0</td><td>[1]</td><td>0</td></tr><tr><td>A20t → 0</td><td>C3</td><td>[b]</td><td>0.0004</td><td>[1]</td><td>0.0002438</td></tr><tr><td>A20t → A20t + A20</td><td>C4</td><td>[b]</td><td>0.5</td><td>[1]</td><td>0.5807</td></tr><tr><td>A20 → 0</td><td>C5</td><td>[b]</td><td>0.0003</td><td>[1]</td><td>0.0003769</td></tr><tr><td>IKKa:IKBα → IKKa + IkBα</td><td>Dα</td><td>[b]</td><td>0.00125</td><td>[2]</td><td>0.002046</td></tr><tr><td>IKKa:IKBβ → IKKa + IkBβ</td><td>Dβ</td><td>[b]</td><td>0.00175</td><td>[2]</td><td>0.0005609</td></tr><tr><td>IKKa:IKBε → IKKa + IkBε</td><td>Dε</td><td>[b]</td><td>0.00175</td><td>[2]</td><td>0.002142</td></tr><tr><td>IKKa:IKBα:NF-κB → IKKa + IkKBα:NF-κB</td><td>Dα</td><td>[b]</td><td>0.00125</td><td>[2]</td><td>0.002046</td></tr><tr><td>IKKa:IKBβ:NF-κB → IKKa + IkBβ:NF-κB</td><td>Dβ</td><td>[b]</td><td>0.00175</td><td>[2]</td><td>0.000561</td></tr><tr><td>IKKa:IKBε:NF-κB → IKKa + IkBε:NF-κB</td><td>Dε</td><td>[b]</td><td>0.00175</td><td>[2]</td><td>0.002142</td></tr><tr><td>IKKa:IKBα:NF-κB → IKKa:IKBα + NF-κB</td><td>Eα</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000144</td></tr><tr><td>IKKa:IKBβ:NF-κB → IKKa:IKBβ + NF-κB</td><td>Eβ</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000124</td></tr><tr><td>IKKa:IKBε:NF-κB → IKKa:IKBε + NF-κB</td><td>Eε</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000064</td></tr><tr><td>IKKa:IKBα + NF-κB → IKKa:IKBα:NF-κB</td><td>Fα</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.3789</td></tr><tr><td>IKKa:IKBβ + NF-κB → IKKa:IKBβ:NF-κB</td><td>Fβ</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.2135</td></tr><tr><td>IKKa:IKBε + NF-κB → IKKa:IKBε:NF-κB</td><td>Fε</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.3528</td></tr><tr><td>IKBα:NF-κB → NF-κB + IkBα</td><td>Gα</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000064</td></tr><tr><td>IKBβ:NF-κB → NF-κB + IkBβ</td><td>Gβ</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000044</td></tr><tr><td>IKBε:NF-κB → NF-κB + IkBε</td><td>Gε</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000069</td></tr><tr><td>IKBα:NF-κBn → NF-κBn + IkBαn</td><td>Gα</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000064</td></tr><tr><td>IKBβ:NF-κBn → NF-κBn + IkBβn</td><td>Gβ</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000044</td></tr><tr><td>IKBε:NF-κBn → NF-κBn + IkBεn</td><td>Gε</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000069</td></tr><tr><td>IKBα + NF-κB → IkBα:NF-κB</td><td>Hα</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.4593</td></tr><tr><td>IKBβ + NF-κB → IkBβ:NF-κB</td><td>Hβ</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.7753</td></tr><tr><td>IKBε + NF-κB → IkBε:NF-κB</td><td>Hε</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.2895</td></tr><tr><td>IKBα + NF-κBn → IKBα:NF-κBn</td><td>Hα</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.4593</td></tr><tr><td>IKBβ + NF-κBn → IKBβ:NF-κBn</td><td>Hβ</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.7753</td></tr><tr><td>IKBε + NF-κBn → IKBε:NF-κBn</td><td>Hε</td><td>[a]</td><td>0.5</td><td>[2]</td><td>0.2895</td></tr><tr><td>NF-κB → NF-κBn</td><td>I1</td><td>[b]</td><td>0.0025</td><td>[1]</td><td>0.003037</td></tr><tr><td>NF-κBn → NF-κB</td><td>K01</td><td>[b]</td><td>0.00005</td><td>[3]</td><td>0.00005537</td></tr><tr><td>IKKn → IKKa</td><td>K1</td><td>[b]</td><td>0.0025</td><td>[1]</td><td>0.003273</td></tr><tr><td>A20 + IKKa → A20 + IKKi</td><td>K2</td><td>[a]</td><td>0.1</td><td>[1]</td><td>0.07075</td></tr><tr><td>IKKa → IKKi</td><td>K3</td><td>[b]</td><td>0.0015</td><td>[1]</td><td>0.00202</td></tr><tr><td>0 → IKKn</td><td>Kprod</td><td>[c]</td><td>0.000025</td><td>[1]</td><td>0.000009752</td></tr><tr><td>IKKn, IKKa, or IKKi → 0</td><td>Kdeg</td><td>[b]</td><td>0.000125</td><td>[1]</td><td>0.0001561</td></tr><tr><td>Volume ratio of cytoplasm to nucleus</td><td>Kv</td><td>1</td><td>5</td><td>[1]</td><td>5</td></tr><tr><td>IKBα:NF-κBn → IkBα:NF-κB</td><td>Lα</td><td>[b]</td><td>0.01</td><td>[1]</td><td>0.013979</td></tr><tr><td>IKBβ:NF-κBn → IkBβ:NF-κB</td><td>Lβ</td><td>[b]</td><td>0.005</td><td>[3]</td><td>0.001567</td></tr><tr><td>IKBε:NF-κBn → IkBε:NF-κB</td><td>Lε</td><td>[b]</td><td>0.005</td><td>[3]</td><td>0.006583</td></tr></tbody></table>
<table><caption>**Table 1** Biochemical reactions and their associated reaction rate constants in the computational model of the NF-κB signaling network (Continued)</caption><tr><td>IκBa:NF-κB → NF-κB</td><td>Mα</td><td>[b]</td><td>0.000025</td><td>[1]</td><td>0.00002837</td></tr><tr><td>IκBβ:NF-κB → NF-κB</td><td>Mβ</td><td>[b]</td><td>0.000025</td><td>[3]</td><td>0.00003609</td></tr><tr><td>IκBε:NF-κB → NF-κB</td><td>Mε</td><td>[b]</td><td>0.000025</td><td>[3]</td><td>0.00000866</td></tr><tr><td>Total NF-κB concentration</td><td>NF-κB</td><td>[d]</td><td>0.06</td><td>[1]</td><td>0.06</td></tr><tr><td>IKKa:IκBa:NF-κB → IKKa + NF-κB</td><td>Pα</td><td>[b]</td><td>0.1</td><td>[1]</td><td>0.12928</td></tr><tr><td>IKKa:IκBβ:NF-κB → IKKa + NF-κB</td><td>Pβ</td><td>[b]</td><td>0.05</td><td>[3]</td><td>0.06454</td></tr><tr><td>IKKa:IκBε:NF-κB → IKKa + NF-κB</td><td>Pε</td><td>[b]</td><td>0.05</td><td>[3]</td><td>0.08434</td></tr><tr><td>IκBan → IκBa</td><td>Qα</td><td>[b]</td><td>0.0005</td><td>[1]</td><td>0.0005123</td></tr><tr><td>IκBβn → IκBβ</td><td>Qβ</td><td>[b]</td><td>0.0005</td><td>[3]</td><td>0.0007398</td></tr><tr><td>IκBεn → IκBε</td><td>Qε</td><td>[b]</td><td>0.0005</td><td>[3]</td><td>0.0002184</td></tr><tr><td>IKKa:IκBa → IKKa</td><td>Rα</td><td>[b]</td><td>0.1</td><td>[1]</td><td>0.123</td></tr><tr><td>IKKa:IκBβ → IKKa</td><td>Rβ</td><td>[b]</td><td>0.1</td><td>[3]</td><td>0.03837</td></tr><tr><td>IKKa:IκBε → IKKa</td><td>Rε</td><td>[b]</td><td>0.1</td><td>[3]</td><td>0.1571</td></tr><tr><td>IκBan:NF-κBn → NF-κBn</td><td>Sα</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.00000037</td></tr><tr><td>IκBβn:NF-κBn → NF-κBn</td><td>Sβ</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.000001131</td></tr><tr><td>IκBεn:NF-κBn → NF-κBn</td><td>Sε</td><td>[b]</td><td>0.000001</td><td>[2]</td><td>0.000001037</td></tr><tr><td>NF-κBn → NF-κBn + IκBat</td><td>Uα</td><td>[b]</td><td>0.0000005</td><td>[1]</td><td>0.000000279</td></tr><tr><td>NF-κBn → NF-κBn + IκBβt</td><td>Uβ</td><td>[b]</td><td>0</td><td>[2]</td><td>0</td></tr><tr><td>NF-κBn → NF-κBn + IκBεt</td><td>Uε</td><td>[b]</td><td>0.00000005</td><td>[3]</td><td>0.000000059</td></tr><tr><td>IκBa → IκBan</td><td>Vα</td><td>[b]</td><td>0.001</td><td>[1]</td><td>0.0009786</td></tr><tr><td>IκBβ → IκBβn</td><td>Vβ</td><td>[b]</td><td>0.001</td><td>[3]</td><td>0.0004871</td></tr><tr><td>IκBε → IκBεn</td><td>Vε</td><td>[b]</td><td>0.001</td><td>[3]</td><td>0.00147</td></tr><tr><td>IκBa, IκBan → 0</td><td>Wα</td><td>[b]</td><td>0.0001</td><td>[1]</td><td>0.000132</td></tr><tr><td>IκBβ, IκBβn → 0</td><td>Wβ</td><td>[b]</td><td>0.0001</td><td>[3]</td><td>0.000133</td></tr><tr><td>IκBε, IκBεn → 0</td><td>Wε</td><td>[b]</td><td>0.0001</td><td>[3]</td><td>0.000042</td></tr><tr><td>IκBat → IκBat + IκBa</td><td>Xα</td><td>[b]</td><td>0.5</td><td>[1]</td><td>0.4552</td></tr><tr><td>IκBβt → IκBat + IκBβ</td><td>Xβ</td><td>[b]</td><td>0.5</td><td>[3]</td><td>0.3828</td></tr><tr><td>IκBεt → IκBat + IκBε</td><td>Xε</td><td>[b]</td><td>0.5</td><td>[3]</td><td>0.3304</td></tr><tr><td>0 → IκBat</td><td>Yα</td><td>[c]</td><td>0.00000005</td><td>[3]</td><td>0.000000084</td></tr><tr><td>0 → IκBβt</td><td>Yβ</td><td>[c]</td><td>0.000000005</td><td>[3]</td><td>0.00000000414</td></tr><tr><td>0 → IκBεt</td><td>Yε</td><td>[c]</td><td>0.000000005</td><td>[3]</td><td>0.00000000508</td></tr><tr><td>IκBat → 0</td><td>Zα</td><td>[b]</td><td>0.0004</td><td>[1]</td><td>0.0003375</td></tr><tr><td>IκBβt → 0</td><td>Zβ</td><td>[b]</td><td>0.0004</td><td>[3]</td><td>0.0002031</td></tr><tr><td>IκBεt → 0</td><td>Zε</td><td>[b]</td><td>0.0004</td><td>[3]</td><td>0.0004742</td></tr></table>

The far left column denotes biochemical reactions. The symbols for reaction rate constants are in the column I. Their units and their published nominal values are in columns II and III respectively. Column IV denotes the reference: [1] for Ref. [21], [2] for Ref. [20], and [3] for the average of values from Refs. [20,21]. Column V lists values used in this paper. The units for [a] are µM⁻¹ s⁻¹, for [b] are s⁻¹, for [c] are µM s⁻¹, and for [d] are µM.

ensemble of single-cell simulations, and compared those distributions for the wild type and various mutants. We made two key findings. First, there is significant overlap between the distributions of the wild type and mutants. This indicates that two individual cells, even if they are genetically different, can respond to the same stimulus in a similar manner. A better way to characterize the differences induced by the differing genetic conditions is to model a large ensemble of cells and compare the full distributions of single-cell responses. Second, for this biochemical pathway, we observed that distributions of

the first Maximum response were the same for any genetic conditions. Similarly, the distributions of the first translocation time responses were the same for the wild type and two of the genetic mutants. This means that some dynamic features are not good indicators of changes in the NF-κB signaling system for genetic comparative studies. The SE approach can be used to screen out bad indicators among the many possible candidates. In the next subsection, we investigate the distributions of dynamic responses for the NF-κB signaling system under two different dosage conditions.
<figure><img src="image_11.png" /><figcaption>**Figure 9** Individual time-series curves (blue lines) and the ensemble average (red) of the key protein concentrations for an ensemble of 1000 replicates of the wild type, stimulated by small dosage (A, C, E, G, and I) or large dosage (B, D, F, H, and J).</figcaption></figure>
## Statistical ensemble analysis of dosage-dependent NF-κB dynamical behavior

### Dosage-dependent dynamical behaviors of individual and ensemble-averaged temporal profiles

We numerically investigated the NF-κB signaling network in response to two different stimulation dosages. As shown in Figure 9, even though a small dosage (TR = 0.01) is 100 times smaller than a large dosage (TR = 1), the small dosage still triggers a substantial amount of NF-κB response from about half the replicates. The other half do not respond at all to the small dosage. In contrast, the large dosage induces strong NF-κB response from all the replicates homogeneously. For example, in Figure 9(A) where the ensemble receives a small dosage, half the temporal profiles of nuclear NF-κB have a single peak and the other half do not. For the half with a peak, the peaks occur after hours of time-delay and there is a large variation in the delays. The steady-state level of nuclear NF-κB concentration is broadly

distributed between zero and 100 nM. In Figure 9(B), the large dosage induces a synchronized appearance of the first peak in all the temporal profiles after a time delay of half an hour. However, even with a large dosage, there is a large cell-to-cell variation both in the amplitude of the first peak and in the timing of successive peaks. In Figures 9(C) and (D), IKK responses to the small dosage stimulation are sharply different from those with large dosage stimulation, i.e., very low levels of IKK versus single-peaked responses with a large amplitude. These IKK profiles are inversely correlated with the profiles of cytoplasmic IκBα. The steady state levels of A20 are 2-3 times higher.

### Dosage-dependent distribution of the dynamic features

The distributions of responses for a thousand-replicate ensemble to large (TR = 1) and small (TR = 0.01) dosage are shown in Figure 10. In Figure 10(A) and (C), both

<figure><img src="image_12.png" /><figcaption>**Figure 10** Distributions of six dynamic features of nuclear NF-κB profiles for an ensemble of 1,000 replicates of the wild type NF-κB signaling system undergoing small (TR = 0.01; red line) or large (TR = 1; black line) dosage stimulations. The six dynamic features are First Maximum in panel (A), First Translocation Time in panel (B), First Period in panel (C), Ratio of First Minimum to First Maximum in panel (D), Ratio of Second Maximum to First Maximum in panel (E), and Ratio of Steady State to First Maximum in panel (F).</figcaption></figure>
the First Maximum and the First Period share similar dosage-dependent behavior: the strength and duration of the response increase with dosage. But, for the First Translocation Time and the Ratio of the First Minimum to the First Maximum metrics in Figure 10(B) and (D), the dosage-dependent behavior is inverted: the larger dosage induces a peak at an earlier time with a smaller First Minimum level. Moreover, the larger dosage makes the distribution more narrowly-distributed. This indicates the larger dosage induces an earlier and spikier response, while the smaller dosage induces more heterogeneous First Maximum and First Minimum levels of nuclear NF-κB concentration. Lastly, both the ratios of the Second Maximum to the First Maximum and of the Steady State to the First Maximum share similar dosage-dependent behavior in Figures 10(E) and 10(F): the smaller dosage induces a distribution at larger values, i.e., closer to one. In other words, when stimulated by the smaller dosage, the levels of the First Maximum, of the subsequent maxima, and of the Steady State are the same, i.e., NF-κB profiles exhibit either a monotonically-increasing pattern or single-peaked pattern with low peak amplitude. In addition, the full half-maximum width of the distribution is unaffected by the dosage amount.

## Dosage-dependent distribution of the dynamic patterns

As shown in Figure 11, when stimulated by a small (TR = 0.01) dosage, 80% of the nuclear NF-κB profiles are damped-oscillatory whereas only 20% of are single-peaked. But, those damped oscillatory responses are similar to a single-peaked response. The distribution induced by the large dosage (TR = 1) corresponds to that of the wild type case in Figure 6. We note that for small dosage stimulation the distribution of the dynamic patterns, the SE average, and the individual profiles of nuclear NF-κB concentration, as shown in Figures 9 and 11, are very similar to those for the mutants responding to large dosage

<figure><img src="image_13.png" /><figcaption>**Figure 11** Distribution of the dynamic patterns of nuclear NF-κB concentration profiles for an ensemble of 1,000 replicates of the wild type NF-κB signaling system undergoing small (TR = 0.01) or large (TR = 1) dosage stimulations. The same four dynamic patterns and coloring scheme are used as in Figure 6.</figcaption></figure>

stimulation with IκBi and IκBε genes double knocked-out, as shown in Figures 5 and 6. We also observed that when the heterogeneity factor χ is increased from χ = 30% to χ = 70%, small dosage stimulation generates more heterogeneous dynamic patterns, i.e. nuclear NF-κB profiles in all the pattern categories.

In this subsection, we analyzed the dynamical response of the cellular replicates under two different stimulant dosage conditions. This yielded distributions of six dynamic features and associated dynamic patterns that are descriptive characterizations of the NF-κB signaling system. Unlike the earlier analysis of differential genetic conditions, the differing stimulus dosages generate non-overlapping distributions and clearly distinctive dynamical behaviors. Some of our predictions, e.g., the distribution of first translocation time in Figure 10(B) and dynamic patterns in Figure 11, are experimentally validated [36].

## Sigmodially shaped SE average of the dose-response curves

We numerically investigated the distribution of dose-response curves from the SE of the NF-κB system. In this analysis we used only 50 replicates of the NF-κB system because of the high computational cost of calculating a single dose-response curve. Each replicate of the NF-κB signaling system is stimulated with a persistent signal for 30 hours, and the average (quasi-steady-state) level of nuclear NF-κB concentration is measured between 20 and 30 hours after stimulation. To check for hysteresis effects, we computed the dose response curve twice, first by increasing the signal strength from TR = 0 to TR = 0.1 in a step-like manner and then by decreasing it from TR = 0.1 to TR = 0. If the forward and backward dose-response curves were significantly different, it could be regarded as a sign of hysteresis. In Figure 12, both forward and backward dose-response curves for each replicate look the

<figure><img src="image_14.png" /></figure>
same, i.e. no hysteresis effects were seen. Figure 12 also
shows both the individual dose-response curves and the
SE average have a sigmoidal shape, which indicates a
switching behavior between on and off states. For a signal
strength greater than the inflection point of a dose-
response curve, the stationary level quickly reaches a pla-
teau whose value is three orders of magnitude higher than
its low stationary level at a signal smaller than the inflec-
tion point. Lastly, the cell-to-cell variability in the statio-
nary nuclear NF-κB level is dramatically larger with a
weak signal than with a strong signal, i.e., the variation is
two orders of magnitude at TR = 10⁻⁴ and one order of
magnitude at TR = 0.1. The variability is a maximum at
the inflection point of the SE dose-response curve, where
it is roughly four orders of magnitude. The cell-to-cell
variability in the location of the inflection point (along
the x-axis) is about one order of magnitude.

Discussion and conclusion

In this paper, we have used a novel statistical ensemble
(SE) method to mimic protein dynamics in a population
of cells influenced by extrinsic noise. For our model of the
NF-κB signaling system, we showed that the SE averages
match population-averaged experimental data. The added
value of the SE method is that it can also produce entire
distributions of response, which can potentially be com-
pared to experimental observations at the single-cell level.
The main predictions enabled by the SE method were as
follows: (a) nuclear NF-κB concentration profiles for single
cells are expected to fall into one of several distinct hetero-
geneous dynamic patterns, (b) larger dosages should induce
more oscillatory dynamic patterns of nuclear NF-κB
response, while smaller dosages should primarily induce
single-peaked patterns, (c) larger (smaller) dosages should
make First translocation times more narrowly-distributed
(broadly-distributed) and shift the peak of its distribution to
earlier (later) times, and (d) the shape of dose-response
curves, both at the single-cell and population level, should
be sigmoidal. After making these predictions computa-
tionally, our experimental colleagues used single-cell
fluorescence imaging to monitor NF-κB nucleo-cytoplasmic
translocation dynamics in lipopolysaccharide-insulted
murine macrophage cells, and found that nuclear GFP-RelA
(a subunit of NF-κB dimers) profiles are very heterogeneous.
They also found NF-κB dynamic responses to be much
more heterogeneous and less oscillatory when the stimulant
dosage was smaller. They also stimulated the murine mac-
rophages with two different dosages (1 nM and 1 µM) of E.
Coli lipopolysaccharide and found two distinctly different
distributions of NF-κB translocation time [36]. Thus two
of our predictions have been verified by our collaborators
who are planning to publish the results elsewhere. We
hope our computational analyses will elicit more single-

cell experimental measurement to verify the predicted
dynamic behaviors.

We wish to emphasize that the novelty of our analysis
is not due to its methodology, but rather the viewpoint
we adopt with regard to computational modeling of cellular
response. Most previous modeling efforts have focused on
bifurcation analysis of the response of a dynamical system
as its input kinetic rate constants are varied. This approach
makes a one-to-one correspondence between a form of
dynamic response and a single set of parameter values.
By contrast, the assumption in the SE approach is that
the dynamics of protein response in individual cells is
intrinsically heterogeneous. We assume a population of
cells (the replicates of the SE) does not occupy a single
point, but rather a volume of points in high-dimensional
parameter space. We choose a hybercube sub-volume of
this space (a midpoint and interval size for each dimension)
and sample from it efficiently to assign kinetic parameters
to each replicate in the SE. In this regard, our SE ap-
proach looks similar to sensitivity analysis by using
Latin Hypercube Sampling method. But, our analysis is
not for sensitivity analysis of cell signaling systems to
perturbation of model parameter values, but for gener-
ation of the heterogeneous responses in single cells. As
explained below, we choose the bounds of this sub-
volume by fitting the resulting SE averages to experi-
mentally observed population-level averages. Once this
is done and the averages match, our assumption is that
we can understand the heterogeneous behavior of the
biochemical network at the single-cell level by analyzing
the wealth of distribution data provided by the SE compu-
tations across its set of replicates.

The sigmoidal shape of the dose-response curve reveals
two important properties of NF-κB signaling: its switching
behavior and its monostability (i.e., no hysteresis). The inflec-
tion points of individual sigmoidal curves can be viewed as
activation thresholds for the NF-κB signaling pathway. As
shown in Figure 12, the NF-κB response is quite small for
signal strength below the threshold, while the response in-
creases dramatically (log scale on y-axis) for signal strengths
just above the threshold. Knowing that some NF-κB target
genes are inflammatory cytokines and that over-expressed
inflammatory response is harmful to the host, we can specu-
late that the NF-κB signaling network employs this sigmoidal
dose-response curve to down-regulate excessive inflamma-
tory responses, i.e., to only turn on if the danger level is sig-
nificantly high, otherwise to shut down. We also note that
the amplitude and timing of the first response peak for in-
flammatory cytokines (such as TNFα) are known to be cri-
tical in mediating timely and effective immune response.
This is motivation for measuring the dosage-dependent tran-
sient dynamic response of NF-κB target genes to investigate
the shape of the dose-response. Lastly, TNFα autocrine sig-
naling forms a positive feedback loop in the NF-κB signaling
network and can induce bistability, which could modify our
results indicating monostability.

Our statistical analysis of protein dynamics depends on
how accurately the computationally generated ensemble of
the NF-κB signaling system represents a true biological
population of individual cells. This question is equivalent to
what is the true distribution of extrinsic noise? I.e., what is
the distribution from which the kinetic input parameters
should be sampled? In this paper, we have chosen a simple
answer to this question by assuming the distribution is uni-
form and bounded. We devised a heuristic fitting algorithm
to find the bounding limits of the uniform distribution for
each kinetic parameter by minimizing the discrepancy be-
tween the SE averages and the population-level experimen-
tal data. This heuristic scheme could be converted to a
more rigorous optimization problem: to find the distribu-
tion of kinetic parameters for a network model which
minimizes the difference between the SE average and the
population-level experimental measurements, while sim-
ultaneously reproducing the range of experimentally-
observed heterogeneous protein dynamics in single cells.

Additional improvements could also be made to the
procedure for sampling from the parameter space. For
example, the sampling could become more biologically
relevant, by accounting for changes in the distribution
of extrinsic noise over time as cells traverse their cell
cycle. We have also assumed no correlation between
pairs of kinetic parameters. In fact some parameters may
be co-dependent because cellular energy resources are
limited: e.g., as one kinetic process is accelerated, others
may be inhibited to balance cellular energy consumption.
All of these computational tasks would be made easier
with additional single-cell experimental data from which
the true distribution of extrinsic noise could be inferred.

Finally, we note that our analysis in this paper was sim-
plified by categorizing the nuclear NF-κB response profiles
into four dynamic patterns. This simplified various statis-
tical analyses and made it easier to characterize changes in
the distribution when genes were knocked-out. Our choice
was based on mathematical characterization of the dy-
namic protein profiles. However, it is possible this neglects
other biologically important details of the nuclear NF-κB
response, e.g. classification by time periodicity or by
steady-state level. Since the choice of categories can affect
subsequent analysis, this is an important factor to consider
when using the SE methodology.

**Methods**

Six dynamic features of nuclear NF-κB profiles

We define six dynamic features to represent the
distinguishing characteristics of temporal profiles of nuclear
NF-κB concentration. The first translocation time is the
time when the first peak occurs; the first period measures
the time between the first two peaks; the first and second

maxima are the amplitudes of the first and second peaks;
the first minimum is the amplitude of the “valley” between
the first two peaks; steady state refers to the asymptotic
amplitude at sufficiently long time. Using the first max-
imum as a reference level, we use scaled ratios, i.e., the first
minimum, the second maximum, and the steady state are
normalized by the first maximum. The distributions of
these dynamical features are presented in Figures 8 and 10.

Generation of the SE of NF-κB signaling network

Each kinetic rate constant listed in Table 1 is randomly
sampled from an interval (x₀(1-χ), x₀(1 + χ)) where χ₀
is the reference value and χ is a heterogeneity factor.
To sample efficiently in the high dimensional space of
dozens of parameters, we use the Latin Hypercube
Sampling methodology discussed below. For this paper,
we used χ = 0.3. To generate a statistical ensemble (SE)
of *N* replicates, we simply generate *N* sets of randomly
sampled kinetic parameters.

Algorithm to fit the SE average to population-level
experimental data

The goal of our fitting algorithm is to determine kinetic
parameters that provide the best match for features of the
SE average to the experimental time-series data. We do not
attempt to fit all of the eighty kinetic parameters. This can
result in “over-fitting”, with too many parameters fit to too
little data. Also, by sensitivity analysis, others and we have
found there are only a handful of kinetic parameters in the
NF-κB signaling network whose variation significantly af-
fects the temporal profile of the nuclear NF-κB concentra-
tion [23,35,40]. Based on our previous studies [35,40], we
choose the two kinetic parameters most highly correlated
with each dynamic feature and varied that set of parameters
in the fitting procedure. The heuristic steps are as follows.
(1) Use an educated guess for initial kinetic parameters and
set the heterogeneity factor to χ = 039. (2) Generate the SE
and resulting protein profiles and calculate the deviation of
the six dynamic features of the SE average from the target
(experimental) dynamic features. (3) Identify the dynamic
feature with largest deviation and modify the two kinetic
parameters associated with it. (4) Repeat steps 1-3 until
the dynamic features are close to the target values. (5)
When a good fit is not achievable decrease χ in a step-like
manner. All of the data in Figures 3, 5, 4, 7 were obtained
through this process.

Numerical simulation of the NF-κB signaling network

A coupled system of ordinary differential equations (ODEs) is derived from the NF-κB signaling network in Figure 2. Using a 4^th^ order Runge-Kutta scheme, we numerically integrate the ODEs, using initial conditions (i.e., the cytoplasmic NF-κB concentration being equal to the total NF-κB concentration and zero concentrations
of all other biochemical species) and kinetic parameters shown in Table 1 and sampled as described below. Before stimulation, the system runs for 33 hours until its constituents reach equilibrium values. Then, we simulate persistent stimulation by turning on the reaction, IKKn → IKKα with a rate TR\*K1 and a non-zero constant value of TR. The ChemCell software package is used to carry out part of numerical simulation [41].

## Latin hypercube sampling (LHS)

LHS is a constrained Monte Carlo sampling scheme. Monte Carlo sampling is a commonly-used approach for assessing the uncertainty of a computational model. By sampling repeatedly from the assumed joint probability function of the input variables, and evaluating the response for each sample, the distribution of responses of the model can be estimated. This approach yields reasonable estimates for the distribution of responses, but a large number of samples is needed if there are many input variables, i.e. a high-dimensional input space. A large sample size can be computationally expensive, which motivates an alternative approach, namely LHS. LHS yields more precise estimates of the response distribution with a smaller number of samples from high-dimensional input spaces [42]. Suppose that the model has K kinetic rate variables and we want N samples, where a sample is a set of K values, one per variable. LHS first selects N different values for each of the K variables, by dividing the range of each variable into N non-overlapping intervals on the basis of equal probability. One value from each interval is selected randomly, in accord with the assumed probability density within the interval. The N values for the first kinetic rate variable are then paired in a random manner (equally likely combinations) with the N values of the second variable. These N pairs are combined in a random manner with the N values of the third variable to form N triplets, and so on, until N K-tuples are formed. Each K-tuple becomes a set of kinetic rate parameters for one replicate within the statistical ensemble of N replicates.

## Competing interests

The authors declare that they have no competing interests.

## Authors' contributions

JJ conceived and designed the numerical experiments and analyzed the simulation results. JJ and SJP developed code specific for these models and performed the numerical experiments. JJ wrote the paper and SJP and JLF edited the paper. All authors read and approved the final manuscript.

## Acknowledgments

This work was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories, a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy, under Contract No. DEAC04-94AL85000. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915.

## Author details

Received: 20 November 2012 Accepted: 7 May 2013

Published: 7 June 2013

## References

1. Thattai M, van Oudenaarden A: Intrinsic noise in gene regulatory networks. Proc Natl Acad Sci 2001, 98:8614-8619.

2. Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Science 2002, 297:1183-1186.

3. Swain PS, Elowitz MB, Siggia ED: Intrinsic and extrinsic contributions to stochasticity in gene expression. Proc Natl Acad Sci 2002, 99:12795-12800.

4. Blake WJ, Karn M, Cantor CR, Collins JJ: Noise in eukaryotic gene expression. Science 2003, 422:633-637.

5. Raser JM, O'Shea EK: Control of stochasticity in eukaryotic gene expression. Science 2004, 304:1811-1814.

6. Rosenfeld N, Young JW, Alon U, Swain SS, Elowitz MB: Gene regulation at the single-cell level. Science 2005, 307:1962-1965.

7. Pedraza JM, van Oudenaarden A: Noise propagation in gene networks. Science 2005, 307:1965-1969.

8. Raj A, Perskin CS, Tranchina D, Vargas DY, Tyagi S: Stochastic mRNA synthesis in mammalian cells. PLoS Biol 2006, 4:e309.

9. Cohen AA, Geva-Zatorsky N, Eden E, Frenkel-Morgenstern M, Issaeva I, Sigal A, Milo R, Cohen-Saidon C, Liron Y, Kam Z, Cohen L, Danon T, Perzov N, Alon U: Dynamic proteomics of individual cancer cells in response to a drug. Science 2008, 322:1511-1516.

10. Spencer SL, Gaudet S, Albeck JG, Burke JM, Sorger PK: Non-genetic origins of cell-to-cell variability in TRAIL-induced apoptosis. Nature 2009, 459:428-432.

11. McAdams HH, Shapiro L: Circuit simulation of genetic networks. Science 1995, 269:650-656.

12. Arkin A, Ross J, McAdams HH: Stochastic kinetic analysis of developmental pathways bifurcation in phage λ-infected E. Coli cells. Genetics 1998, 149:1633-1648.

13. Blake WJ, Balazsi G, Kohanski MA, Isaacs FJ, Murphy KF, Kuang Y, Cantor CR, Walt DR, Collins JJ: Phenotypic consequences of promoter-mediated transcriptional noise. Mol Cell 2006, 24:853-865.

14. Paulsson J: Summing up the noise in gene networks. Nature 2004, 427:415-418.

15. Shahrezaei V, Ollivier JF, Swain PS: Colored extrinsic fluctuations and stochastic gene expression. Mol Syst Biol 2008, 4:196.

16. Brown KS, Hill CC, Calero GA, Myers CR, Lee KH, Sethna JP, Cerione RA: The statistical mechanics of complex signaling networks: nerve growth factor signaling. Phys Biol 2004, 1:184-195.

17. Huang K: Statistical Mechanics. secondth edition. New York: Wiley; 1987.

18. Verma IM, Stevenson J: IkB kinase: beginning, not the end. Proc Natl Acad Sci U S A 1997, 94:11758-11760.

19. Li Q, Verma IM: NF-κB regulation in the immune system. Nat Rev Immunol 2002, 2:725.

20. Hoffmann A, Baltimore D: Circuitry of nuclear factor κB signaling. Immunol Rev 2006, 210:171-186.

21. Hoffmann A, Levchenko A, Scott ML, Baltimore D: The IkB-NF-κB signaling module: temporal control and selective gene activation. Science 2002, 298:1241-1245.

22. Lee EG, Boone DL, Chai S, Libby SL, Chien M, Lodolce JP, Ma A: Failure to regulate TNF-induced NF-κB and cell death responses in A20-deficient mice. Science 2000, 289:2350-2354.

23. Lipniacki T, Paszek P, Brasier AR, Luxon BA, Kimmel M: Mathematical model of NF-κB regulatory module. J Theor Biol 2004, 228:195-215.

24. Lipniacki T, Paszek P, Brasier AR, Luxon BA, Kimmel M: Stochastic regulation in early immune response. Biophys J 2006, 90:725-742.

25. Nelson DE, Ihekwaba AEC, Elliott 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 2004, 306:704-708.

26. Nelson DE, Horton CA, See V, Johnson JR, Nelson G, Spiller DG: Response to comments on "oscillations in NF-κB signaling control the dynamics of gene expression". Science 2005, 308:52b.

<footnote n="1">¹Department of Physics and Astronomy, University of Tennessee, Knoxville 37996, USA. ²Scalable Algorithms Department, Sandia National Laboratories,</footnote>

<footnote>Albuquerque, NM 87185, USA. ³Department of Biology, Evry University, Evry Cedex, France.</footnote>
27. Barken D, Wang CJ, Kearns J, Cheong R, Hoffmann A, Levchenko A: Comment on "oscillations in NF-κB signaling control the dynamics of gene expression". *Science* 2005, **308**:52a.

28. Lee TK, Denny EM, Sanghvi JC, Gaston JE, Maynard ND, Hughey JJ, Covert MW: A noisy paracrine signal determinines the cellular NF-κB response to lipopolysaccharide. *Sci Signal* 2009, **2**(93):ra65.

29. Ashall L, Horton CA, Nelson DE, Paszek P, Harper CV, Sillitoe K, Ryan S, Spiller DG, Untt JF, Broomhead DS, Kell DB, Rand DA, Sée V, White MR: Pulsatile stimulation determines timing and specificity of NF-κB-dependent transcription. *Science* 2009, **324**:242-246.

30. Bartfeld S, Hess S, Bauer B, Machuy N, Ogilvie LA, Schuchhardt J, Meyer TF: High-throughput and single-cell imaging of NF-κB oscillations using monoclonal cell lines. *BMC Cell Biol* 2010, **11**:21.

31. Turner DA, Paszek P, Woodcock DJ, Nelson DE, Horton CA, Wang Y, Spiller DG, Rand DA, White MRH, Harper CV: Physiological levels of TNF stimulation induce stochastic dynamics of NF-κB responses in single living cells. *J Cell Sci* 2010, **123**:2834-2843.

32. Tay S, Hughey JJ, Lee TK, Lipniacki T, Quake SR, Covert MW: Single-cell NF-κB dynamics reveal digital activation and analogue information processing. *Nature* 2010, **466**:267-271.

33. Paszek P, Ryan S, Ashall L, Sillitoe K, Harper CV, Spiller DG, Rand DA, White MR: Population robustness arising from cellular heterogeneity. *Proc Natl Acad Sci* 2010, **107**:11644-11649.

34. Hayot F, Jayaprakash C: NF-κB oscillations and cell-to-cell variability. *J Theor Biol* 2006, **240**:583-591.

35. Ihekwa AEC, Broomhead DS, Grimley RL, Benson N, Kell DB: Sensitivity analysis of parameters controlling oscillatory signaling in the NF-κB pathway: the roles of IKK and IκBα. *Syst Biol* 2004, **1**:93-103.

36. James CD, Moorman MW, Carson BD, Branda CS, Lantz JW, Manginell RP, Martino A, Singh AK: Nuclear translocation kinetics of NF-κB in macrophages challenged with pathogens in a microfluidic platform. *Biomed Microdevices* 2009, **11**:693-700.

37. Tian B, Nowak DE, Brasier AR: A TNF-induced gene expression program under oscillatory NF-κB control. *BMC Genomics* 2005, **6**:137-155.

38. Werner SL, Barken D, Hoffmann A: Stimulus specificity of gene expression programs determined by temporal control of IKK activity. *Science* 2005, **309**:1857-1861.

39. Covert MW, Leung TH, Gaston JE, Baltimore D: Achieving stability of lipopolysaccharide-induced NF-κB activation. *Science* 2005, **309**:1854-1857.

40. Joo J, Plimpton S, Martin S, Swiler L, Faulon JL: Sensitivity analysis of a computational model of the IKK-NF-kappaB-IkappaB-alpha-A20 signal transduction network. *Ann NY Acad Sci* 2007, **1115**:221-239.

41. Plimpton SJ, Slepoy A: Microbial cell modeling via reacting diffusing particles. *J Phys: Conference Series* 2005, **16**:305-309.

42. Swiler LP, Wyss GD: A user's guide to Sandia's Latin Hypercube sampling software: LHS UNIX library/standalone version, SAND Report. 2004:2004-2439.

doi:10.1186/1752-0509-7-45

Cite this article as: Joo et al.; Statistical ensemble analysis for simulating extrinsic noise-driven response in NF-κB signaling networks. BMC Systems Biology 2013 7:45.

Submit your next manuscript to BioMed Central and take full advantage of:

* Convenient online submission

* Thorough peer review

* No space constraints or color figure charges

* Immediate publication on acceptance

* Inclusion in PubMed, CAS, Scopus and Google Scholar

* Research which is freely available for redistribution

Submit your manuscript at
www.biomedcentral.com/submit

( ) BioMed Central