Sources of Cell-to-cell Variability in Canonical Nuclear Factor-κB (NF-κB) Signaling Pathway Inferred from Single Cell Dynamic Images*

The canonical nuclear factor-κB (NF-κB) signaling pathway controls a gene network important in the cellular inflammatory response. Upon activation, NF-κB/RelA is released from cytoplasmic inhibitors, from where it translocates into the nucleus, subsequently activating negative feedback loops producing either monophasic or damped oscillatory nucleo-cytoplasmic dynamics. Although the population behavior of the NF-κB pathway has been extensively modeled, the sources of cell-to-cell variability are not well understood. We describe an integrated experimental-computational analysis of NF-κB/RelA translocation in a validated cell model exhibiting monophasic dynamics. Quantitative measures of cellular geometry and total cytoplasmic concentration and translocated RelA amounts were used as priors in Bayesian inference to estimate biophysically realistic parameter values based on dynamic live cell imaging studies of enhanced GFP-tagged RelA in stable transfectants. Bayesian inference was performed on multiple cells simultaneously, assuming identical reaction rate parameters, whereas cellular geometry and initial and total NF-κB concentration-related parameters were cell-specific. A subpopulation of cells exhibiting distinct kinetic profiles was identified that corresponded to differences in the IκBα translation rate. We conclude that cellular geometry, initial and total NF-κB concentration, IκBα translation, and IκBα degradation rates account for distinct cell-to-cell differences in canonical NF-κB translocation dynamics.

sis. In addition, it was observed that NF-B also potently activates expression of IB␣ and TNFAIP3/A20, two proteins mediating discrete steps of the negative feedback regulation of the NF-B signaling pathway. IB␣ expression recaptures NF-B/RelA back into the cytoplasm, regenerating the cytoplasmic IB␣ pool (14,15), whereas the deubiquitinase TNFAIP3/A20 binds and inactivates IKK (16), terminating IB␣ proteolysis. Thus, NF-B regulates multiple inhibitory feedback loops that ensure only transient activation of inflammation (17).
The understanding of dynamic regulation of the canonical signaling pathway has been enhanced by theoretical approaches. Deterministic ordinary differential equation (ODE) models have been developed using population-derived data to examine the specific functional roles of IB␣, -␤, and -⑀ isoforms in regulating the temporal control of NF-B (18) and the role of the NF-B-TNFAIP3/A20 feedback loop (19) (reviewed in Ref. 20). New advances in dynamic microscopy have enabled single cell imaging experiments (21)(22)(23). These have supplemented population-based experiments to indicate that the dynamic behavior of NF-B nuclear translocation can vary significantly between cells. In response to tonic TNF␣ stimulation, an initial synchronous pulse of NF-B translocation is followed either by nuclear persistence (monophasic activation) or in a minority of cells by dampened oscillatory profiles (22,24). Individual cells producing asynchronous oscillatory responses produce population robustness, a phenomenon that is predicted to reduce tissue sensitivity to paracrine stimuli (26).
Inspection of the dynamic imaging studies indicates that cell populations exhibit heterogeneity in the kinetics of the initial NF-B response profile. The sources producing cell-to-cell variation are largely unknown. Sensitivity analysis in simulations has indicated that cellular geometry, represented by the ratio of cytoplasmic to nuclear volume, a parameter termed k v , has a significant influence on the initial spike of NF-B because these molecules translocate into the nucleus. Using a clustering approach on trajectory simulations sampled over the parameter space, six different dynamic types of NF-B translocation behavior were identified (27). However, the experimental measurement of cellular geometry in single cell measurements has not yet been accomplished. Others have proposed stochastic responses at the level of gene expression, although these responses largely influence the oscillatory modes (28,29).
Despite theoretical advances in NF-B pathway modeling, individual reaction rates and initial parameters are underdetermined. To extend our understanding of the pathway, we examined dynamic regulation of the NF-B pathway in an experimentally validated model of airway epithelial cells stably expressing an enhanced green fluorescent protein (EGFP)tagged NF-B/RelA molecule. Quantitative measurements of NF-B translocation surprisingly indicate that only a small fraction (Ͻ20%) of cytoplasmic RelA is translocation-competent, being distributed in IB␣, -␤, and -⑀ isoform-containing complexes. Single cell dynamic imaging revealed that RelA translocation dynamic was monophasic in these cells. Cellular geometry was experimentally measured using axial sectioning. Incorporation of the cellular geometry and fractional NF-B translocation into the priors for the Bayesian inference (30) aided the determination of joint probability distribution for the relevant parameters in the translocation model. Our analysis suggests a subpopulation of cells with significantly faster NF-B decay profiles due to differences in the IB␣ translation and degradation dynamics.

EXPERIMENTAL PROCEDURES
Materials-Recombinant TNF␣ (Pepro Tech) was added directly to the culture medium at a final concentration of 30 ng/ml.
Cell Lines and Extracts-The pCX4-pur-EGFP-RelA-expressing plasmid was constructed by cloning EGFP-RelA cDNA into pCX4-pur vector (31). The EGFP-RelA cDNA was isolated from pcDNA-EGFP-RelA (32) by AgeI and XbaI restriction endonuclease digestion and ligated into pCX4-pur vector. Bosc 23 cells were cotransfected with pCX4-pur-EGFP-RelA and the amphotrophic packaging plasmid pCL-10A1 using Lipofectamine 2000 (Invitrogen). After 3 days, retrovirus in the conditioned medium was collected, filtered, and used to infect A549 cells at 30 -50% confluence in the presence of 8 g/ml Polybrene. After 2 days of infection, the culture medium was changed to contain 2 g/ml puromycin for selection. Puromycin-resistant clones were isolated and transferred to 6-well cell culture plates. Only strongly green fluorescent clones were selected further. A549 EGFP-RelA cells were grown in a medium containing 90% Dulbecco's modified Eagle's medium, 10% heat-inactivated fetal bovine serum, 4 mM L-glutamine, 0.1 mM non-essential amino acids, 1 mM sodium pyruvate, 100 mg/ml G418, 100 units/ml penicillin G sodium, and 100 mg/ml streptomycin sulfate in a humidified atmosphere of 5% CO 2 . Cytoplasmic and sucrose cushion nuclear extracts were prepared as described (11,33). Protein concentration was measured using Coomassie Blue staining using bovine serum albumin as a standard (Bio-Rad) or by amino acid analysis for absolute protein quantification.
For confocal fluorescence microscopy, 25-mm round microscope cover glasses (Fisher) were precoated with collagen solution (Roche Applied Science) in a culture plate. A549 EGFP-RelA cells were plated 1 day prior to the experiment. For imaging, the cover glass was mounted into a microscopic incubation chamber where cells were maintained at 37°C with humidified 5% CO 2 .
Time Lapse Microscopy-Dynamic live cell imaging of EGFP-RelA transfected A549 cells was performed using a Zeiss LSM-510 META confocal microscope with a 63ϫ, 1.4 numerical aperture oil immersion objective with an incubation system where cells were stably maintained at 37°C with humidified 5% CO 2 . Samples were excited using the 488-nm line at 1% of transmission, and emission was collected with a low pass 505-nm filter. The time lapse images for each sample were acquired at 6-min intervals. Focus was automatically corrected before each time point by a customized Zeiss Multitime autofocus macro.
Image Quantification and Analysis-Cell images were segmented, quantified, and tracked by a custom written pipeline program using CellTracker (34). In each time lapse field, cells were excluded from analysis that showed overlap with other cells, divided during the time course, or moved out of view. Images were smoothed using a median filter of size six and normalized using local average and its S.D. value. Both cell and nuclear boundaries were then identified first automatically by intensity thresholding followed by manual correction. The same cells were tracked from one frame to the next frame through the time course. Data from 20 single cells were collected for each time course profile. Where indicated, the nuclear EGFP-RelA fluorescence signal was divided by the total cellular EGFP-RelA signal in the measurement volume to obtain the fraction of the NF-B/RelA molecules that translocated.
Sequential Axial Image Collection and Analysis-The cytoplasmic/nuclear volume ratio (k v ) was measured from sequential axial images obtained using a Zeiss LSM-510 META confocal microscope with a 63ϫ, 1.4 numerical aperture oil immersion objective. For each of the TNF␣-stimulated and non-stimulated fixed cells, images were collected using two lines of excitation, 488 and 633 nm (both at 5% of transmission), and two different channels of emission. The intensity of emission for EGFP-RelA was obtained by using a 505-530-nm band pass filter (Ch2-T1), and the intensity of fluorescence of To-Pro3 was obtained by using a 650-nm low pass filter (Ch3-T2), respectively. All images were collected using 4-frame Kalman averaging with a pixel time of 2.51 s and a pixel size of 200 ϫ 200 ϫ 300 nm, hence z-steps of 0.3 m at an optical slice of 0.7 m. Scan zoom was 1.4, and frame size was 512 ϫ 512. Images were imported to Metamorph software (Molecular Devices, Downingtown, PA) for processing. Initially, background was subtracted from all images based on the signal of non-cellular regions. Regions were traced around the images collected on Ch3-T2 (nuclei) and on Ch2-T1 (cell body). Areas of the desired regions were then calculated, and the volume for that particular slice was obtained by multiplying region areas by 0.3 m as the z-dimension. The total volumes of the cell compartments (cells and nuclei) were then computed by summing the volumes for the individual slices.
Western Immunoblots-Cytoplasmic or nuclear extracts from a constant number of cells were boiled in Laemmli buffer, fractionated on 10% SDS-PAGE, and electrotransferred to polyvinylidene difluoride membranes (Millipore, Bedford, MA). Membranes were blocked in 5% milk, Tris-buffered saline, 0.1% Tween (TBS-T) for 1 h and immunoblotted with the indicated primary Ab for 1 h at 4°C (Santa Cruz Biotechnology, Inc.). Membranes were washed four times in TBS-T. Immune complexes were detected by infrared dye (IRD)-labeled goat anti-rabbit Ab and imaged in an Odyssey infrared imager (Li-Cor, Lincoln, NE). Standard curves were constructed for recombinant RelA to estimate RelA abundance. Homogeneous recombinant GST-RelA(488 -551) was quantified by amino acid analysis. Serial dilutions of the purified GST-RelA(488 -551) protein were assayed to generate standard curves for each immunoblot.
Immunoprecipitations-Cytoplasmic extracts were subjected to depleting immunoprecipitation using 4 g of rabbit anti-IB␣, -␤, or -␥ Abs, and immune complexes were captured on protein A-agarose beads (40-l slurry, Sigma-Aldrich). This amount was determined in separate experiments to deplete Ͼ90% of the target protein. The abundance of RelA with each IB isoform was then determined in Western immunoblots using anti-RelA antibody to quantify RelA associated with each IB isoform. Depletion of the respective IB isoform was demonstrated by Western immunoblot of the flow-through of the immunoprecipitate. For the sequential immunoprecipitation experiment, 4 g of rabbit anti-IB␣ Ab was incubated with 2 mg of cytoplasmic lysate. The immune complexes were then removed using 40 ml of Protein A magnetic beads (Dynal, Invitrogen). The supernatant was then sampled for Western immu-noblot for RelA, and the remainder was incubated with 4 mg of anti-IB␤ Ab. The process was repeated for removal of NF-B1 and -2, as indicated.
Mathematical Modeling Reactions and Parameters-The model involves two-compartment kinetics of the NF-B signaling pathway: IKK and NF-B as activator proteins and A20 and IB␣ as inhibitor proteins (Fig. 3). Because all proteins in the system are present in large numbers, as reported under "Results," the effects of stochasticity in individual molecular reactions on the time evolution of each cell are expected to be minor. Therefore, an ODE model is chosen to describe the NF-B translocation kinetics based on Ref. 19. The mathematical representation consists of a system of 14 ODEs and 27 reactions that account for complex formation of NF-B⅐IB␣ and TNFAIP3/A20-IKK as well as IB␣ degradation and RelA and IB␣ nuclear and cytoplasmic transport (see the supplemental material for a listing of the equations and the associated nominal rate constants). The IKK complex and the dimer NF-B are considered as single proteins; hence, the dynamics leading to formation of these complexes will be disregarded. It is assumed that the cytoplasmic complex IKK may exist in one of three forms: neutral (IKKn), active (IKKa), or inactive (IKKi). We also assume that each form of IKK undergoes degradation with the same degradation rate (k deg ); therefore, the total amount of IKK (i.e. IKKn ϩ IKKa ϩ IKKi) remains constant at equilibrium. The IKKa can form transient complexes with IB␣ proteins or IB␣⅐NF-B complexes. Based on Western immunoblot analysis, we have introduced another form of NF-B, termed "reservoir NF-B," which does not translocate into the nucleus. A small percentage of total NF-B that is translocation-competent is termed "active NF-B." The exchange reactions between reservoir NF-B and active NF-B were assumed to be very slow compared with the translocation dynamics, and their rates were therefore set to zero in the current analysis. The mRNA synthesis rate is the sum of a steady term and of an inducible term and is the same for both the A20 and IB␣ proteins. The ODE model is presented in the supplemental material.
Sensitivity Analysis Methodology-Sensitivity analysis was used to rank model variables that were most important in controlling the pathway dynamics. We investigated the sensitivity of 30 parameters using a polynomial chaos representation. The 30 parameters include all 27 reaction rates, initial NF-B concentration (NF-B bound to IB␣), total NF-B concentration, and k v ratio (supplemental material) (35).
Bayesian Inference-Bayesian framework leads to probabilistic estimates of the model parameters in sparse and noisy data sets. The approach relies on Bayes' formula, m͑t͒ ϭ I n ͑t͒ I n ͑t͒ ϩ I c ͑t͒ ϭ N n ͑t͒ N n ͑t͒ ϩ N c ͑t͒ Data D-The measured absolute intensities in nucleus (I n ) and cytoplasm (I c ) are proportional to the number of NF-B molecules in the nucleus (N n ) and cytoplasm (N c ), respectively (i.e. N n ϭ ␣I n and N c ϭ ␣I c ). However, because ␣ is unknown, we consider only the ratio of nuclear intensities to the total ones and relate them to the model variables (i.e. concentrations) in the following way (including time dependence), where NF-B n T (t) and NF-B c T (t) are the total nuclear and total cytoplasmic concentrations for NF-B in any form. To be precise, m(t) is the model output for the nuclear intensity fraction, whereas the corresponding data set D ϭ {d(t i )} of measurements taken at different time points t i , i ϭ 1,2, . . . N forms the data that are used in Bayes' formula.
Model M ϭ { 3 c}ʚ{Rate Constants, Initial Conditions, k v , T c }-The model in this context is a vector of some of the ODE model parameters (i.e. rate constants, initial conditions, the volume ratio k v , and the total concentration, Prior P(M)-In the absence of any specific shape information, we assume a uniformly distributed prior for the logarithms of the rate constants across a reasonable physical range for each of them. The prior information for the total concentration T c was obtained experimentally and assumed to be a log normal distribution. We assume a log normal prior distribution for the volume ratios k v as well, with parameters that are found by fitting the "area" data obtained from z-stack imaging data.
Evidence P(D)-The evidence merely plays a role of normalizing constant and does not affect the posterior estimation, as will become clear below.
Likelihood P(D͉M)-The likelihood construction is viewed as a function of the model M and is equal to the probability of having a particular data set given the model. This conveniently incorporates errors produced by measurement as well as errors associated with the fact that the model is a simplification of the real network. This discrepancy is assumed to be gaussian-distributed and independent for each time snapshot. The underlying gaussian error model is assumed to have zero mean and an S.D. value that is linearly related to the measured data, where i 2 ϭ s p ϩ s q d(t). The linear dependence coefficients, s p and s q , are hyperparameters. In the single cell inference case, the hyperparameters are inferred along with the model parameters, whereas in the multicell inference, we have fixed them (to keep the dimensionality of the Markov chain Monte Carlo (MCMC) algorithm from being too high) to values that are within the range of the posterior estimates found from single cell inference. Namely, these parameters are set to s p ϭ 0.01 and s q ϭ 0.03. The assumed discrepancy model (Equation 4) then leads to the likelihood formula, where N is the number of data points. The product is a direct consequence of the independent measurement assumption. Furthermore, in order to sample from the posterior distribution via Bayes' formula, an MCMC technique has been used, with the original Metropolis algorithm (36) enhanced by adaptive proposal distributions (37) (see the supplemental material). Also, to ensure positive values for all of the parameters, MCMC was performed on the logarithms of the values of the parameters.
Single Cell Inference-Bayesian inference was performed for each cell based on respective cell data. In the single cell analysis, the dimensionality of the MCMC chain is eight parameters (i.e. three reaction kinetic rate constants (c 4a , c 6a , and i 1 ) as well as the volume ratio k v , initial NF-B⅐IB␣ concentration, and total NF-B concentration T c , along with hyperparameters s p and s q , were inferred). We have obtained a reasonable chain mixing and acceptance rate of ϳ50% for N MCMC ϭ 80,000 steps by using an initial gaussian proposal distribution of width 0.05 (i.e. C 0 ϭ 0.05I d ; see supplemental material) for each parameter during the first n o ϭ 5,000 steps and adaptive MCMC covariance scaling factor ␥ ϭ 0.05 for the remaining 75,000 steps.
Simultaneous Inference of Multiple Cells-For all three multicell analyses (i.e. 20-, 16-, and 4-cell analyses), we allowed the reaction kinetic rate constants (c 4a , c 6a , and i 1 ) to be common among cells during the inference, whereas other parameters (k v , initial NF-B, and T c ) were made specific to each cell and were inferred. The rationale behind this strategy is that the probability of a cell property to vary from other cells is higher for its cell volume, its compartmental volume ratios, and its NF-B concentration rather than its physiological reaction rates. In order to avoid a very high dimensional MCMC, the noise model parameters were set to s p ϭ 0.01 and s q ϭ 0.03 for all cells. In 20-cell inference, for instance, this resulted in a 63-dimensional space for the inference problem. The rest of the MCMC parameters were kept the same as in the single cell analysis.

RESULTS
Validation of the Cellular Model-To measure dynamic changes in RelA translocation at the single cell level, we constructed human A549 alveolar epithelial cells stably expressing an EGFP-tagged NF-B/RelA. Because ectopic EGFP-RelA expression may influence the dynamics of the system response, we characterized the cells to determine the relative levels of ectopic to endogenous RelA and to assess translocation dynamics. Western immunoblots were performed of control and stimulated cytoplasmic and sucrose cushion-purified nuclear extracts (these fractionation methods have been extensively validated, with the nuclear fraction free of cytoplasmic contamination, measured by ␤-tubulin staining, and the cytoplasmic fraction free of nuclear contamination, measured by ␤-laminin staining (33, 38)). Here, in unstimulated cell cytoplasm, the 96-kDa EGFP-RelA band could be identified at a ϳ1:1 ratio with that of 65-kDa endogenous RelA in A549 wild type (WT) cells (Fig. 1A). Upon TNF␣ stimulation, both EGFP-RelA and endogenous RelA translocated into the nucleus. Interestingly, TNF␣ stimulation did not result in significant cytoplasmic depletion of RelA in either wild type or EGFP-RelA-expressing cells (Fig. 1A). This behavior is similar to other cell types where the fraction of cytoplasmic protein translocating into the nucleus was analyzed (15).
To confirm these data, comparative confocal imaging experiments were conducted. In unstimulated cells, EGFP-RelA was primarily distributed in the cytoplasm. In response to TNF␣ stimulation, EGFP-RelA translocated into the nucleus (exclusive of nucleoli) with kinetics similar to that observed for the endogenous RelA protein in A549 WT cells (Fig. 1B). Additionally, we noted that although the nuclear staining was increased in response to stimulation, detectable amounts of RelA remained cytoplasmic for both EGFP-RelA and endogenous RelA in A549 WT cells (Fig. 1B).
We and others have shown that NF-B/RelA controls the expression of a network of genes whose expression occurs in temporally distinct waves, with some genes activated rapidly (IB␣, Gro␤, TNFAIP3/A20, IL-8, and IL-6) and others being activated slowly (IB⑀, NF-B2, and Naf-1) in response to TNF␣ stimulation (11,13). To further characterize our cellular model, the profiles of representative members of the NF-Bdependent gene network were analyzed in EGFP-RelA and A549 WT cell response to tonic TNF␣ stimulation. Members of the "early" NF-B-dependent subnetwork were strongly induced by TNF␣, peaking within 1 h of stimulation in a manner that was identical for both the EGFP-RelA and the A549 WT cells (Fig. 1C). Similarly, the "late" subnetwork was activated after 6 -9 h, with kinetics and magnitude that were indistinguishable between EGFP-RelA and A549 WT cells (Fig. 1C). These data indicate that the EGFP-RelA is expressed at nearly endogenous levels and does not detectably affect the kinetics of RelA translocation or expression of downstream NF-B-dependent genes.
Identification of a Cytoplasmic NF-B Reservoir-The qualitative Western immunoblots and the confocal microscopy suggested that RelA remained cytoplasmic, even during peak NF-B translocation in response to a saturating TNF␣ dose (30 ng/ml). To confirm this observation, we fractionated resting and TNF␣-stimulated cells into cytoplasmic and nuclear extracts and measured changes in RelA abundance using quantitative Western immunoblot, where serial dilutions of purified recombinant RelA were co-fractionated for use as a standard curve ( Fig. 2A). From this analysis, we estimated that A549 WT cells express ϳ288,000 molecules of endogenous RelA per cell, and nearly 23,000 translocate, a number constituting ϳ10% of the total cytoplasmic molecules available. These data indicate that a significant fraction of RelA remains cytoplasmic, a finding confirmed using an independent aptamer pull-down technique to quantify activated RelA in these same cells (39).
To determine whether this IB␣/␤/⑀-associated RelA represented the majority of the cytoplasmic RelA pool, we conducted a depleting immunoprecipitation experiment. Cytoplasmic extracts from unstimulated A549 cells were incubated individually with IB␣, -␤, -␥, -⑀, or -or Bcl-3 Abs, and RelA-associated complexes were removed by binding protein A-agarose beads. The remaining, non-IB isoform-complexed RelA was measured in the supernatant by quantitative Western blot (Fig.  2C). Although IB␣ was completely depleted from the cytoplasmic extract (Fig. 2C, middle), the majority of RelA remained in the supernatant, again confirming that a minority of the total cytoplasmic RelA is complexed with IB␣, -␤, and -⑀.
Finally, a sequential immunodepletion experiment was conducted to study the association of RelA with NF-B1/p50, NF-B2/p52, and IB(␣ and ␤) proteins (a schematic representation is shown in Fig. 2D). Quantitation of the images indicates that 48% of RelA is associated with NF-B1/p50 and 15% is associated with NF-B2/p52 (with the remainder of RelA associating with IB isoforms; Fig. 2E). Because NF-B1 and -2 are not proteolyzed by the TNF-induced IKK␤ pathway (12,40), these complexes largely account for the inert cytoplasmic sequestration of RelA. Together, we interpret these data to indicate that a significant fraction of RelA in A549 cells is not associated with IB isoforms and therefore is not translocationcompetent in response to TNF␣.
As an independent estimation of the fraction of NF-B translocation, we performed serial axial images of control and TNF␣-stimulated EGFP-RelA-expressing cells. Integrated fluorescence intensity was summed for the entire cell, and the fraction that was in the nuclear compartment was expressed as a fraction of the total fluorescence intensity. These data also indicated that 8 -14% of the total EGFP-RelA signal translocated into the nucleus after stimulation (Fig.  2F), a finding that confirmed the cellular fractionation experiments ( Fig. 2A).
Measurement of Cytoplasmic/Nuclear Volume Ratio-In mathematical representations of the NF-B pathway, the cellular geometry parameter k v , the ratio of cytoplasmic to nuclear volume, plays an important role in pathway dynamics, determining the nuclear concentration of translocated RelA (19,27). Because this parameter has not been experimentally determined, we performed volume measurements using confocal microscopy. In these measurements, control and TNF-stimulated EGFP-RelA cells were fixed and imaged. Multiple 0.7-mthick sequential images were taken along the vertical (z) dimen-sion with an optical slice separation of 0.3 m. These slices were assembled, and k v values were computed. These data show that the k v values vary from cell to cell, ranging from 2 to 4.5 with a mean of 3. The fitted lognormal prior distribution with this experimental data with 90% confidence interval is shown in Fig.  S1. Moreover, k v is not significantly affected by TNF␣ treatment (see Table 1). Model Formulation-Both confocal and Western blot experiments indicated the existence of a reservoir of non-IB-associated RelA that does not translocate to the nucleus. Based on these convergent experimental results, we modified a two-level negative feedback deterministic model to include "active" and "inactive" NF-B forms (Fig. 3). The inactive form is assumed to be sequestered in a reservoir, with NF-B1 or -2, and is assumed to be in slow exchange with the active translocation-competent RelA, bound by IB complexes. The rate of exchange between these two NF-B reservoirs are assumed to be very small on the time scales of the TNF␣ stimulation experiments and were therefore set to zero. The inactive reservoir and active NF-B together account for the total cellular concentration of NF-B measured by the dynamic imaging experiments.
EGFP-RelA Cells Show Monophasic Translocation Profiles-EGFP-RelA transfectants were imaged using dynamic confocal microscopy. Sequential images were first obtained to establish basal conditions. The cells were then stimulated with a saturating dose of TNF␣ (30 ng/ml), a dose that produces maximal effects on RelA translocation and target gene expression. Serial confocal images were taken every 6 min, and EGFP fluorescence was quantified and expressed as a fractional translocation. Inspection of the time series for images up to 6 h shows that the RelA translocation profile is monophasic (not shown), a finding consistent with the nuclear translocation of endogenous RelA in our previous studies (11). The profile for the first 72 min of stimulation is shown in Fig. 4.
Further analysis of the RelA translocation response clearly indicated that the time until peak nuclear translocation was highly variable, with 25% cells showing a peak translocation between 41 and 50 min after stimulation. The basal and peak NF-B values differed between cells. In addition, we noted distinct profiles of decay of the monophasic peak, examined below. Together these data indicate the presence of cellular variability in the canonical RelA translocation profile in its non-oscillatory activation mode.
Refinement of Initial Concentrations of NF-B as the IB␣⅐NF-B Complex-Initial ODE simulations using previous model formulations assuming that the total NF-B was translocation-competent did not agree with the basal levels of nuclear NF-B for the cells under study. These higher initial  , c 1 and c 1a ). These two proteins act as negative feedback loops to turn off the NF-B signaling (dotted lines and translation reactions, c 4 and c 4a ). Experimentally shown in this study, a large amount of NF-B does not translocate to the nucleus upon stimulation and acts as a reservoir of NF-B. The exchange between the reservoir and translocable NF-B is shown with dashed lines; however, their exact mechanism and rate kinetics (P1, P2) are yet to be determined and are assumed to be relatively slow relative to the overall NF-B translocation kinetics. The total concentration of NF-B is determined by quantitative Western blots in this study to be ϳ0.12 M. The initial translocation-competent fraction of NF-B is represented by NF-B⅐IB␣ complex. See supplemental material for specific parameter values. concentrations resulted in oscillatory behavior more frequently than did low initial concentrations for the same cell geometry (data not shown). This observation suggested an overestimation of the initial NF-B concentration in the translocationcompetent form, represented by the IB␣⅐NF-B complex. Therefore, this initial NF-B concentration was inferred in our modified model along with the total NF-B concentration (including the reservoir NF-B). Prior information on these concentrations was obtained from Fig. 2A (also see lognormal prior distribution, Fig. S2), where the average total NF-B concentration in the cell population was calculated as 0.12 M and the active NF-B in the form of IB␣⅐NF-B complex was estimated to be 0.035 M. We note that this concentration is nearly half the value of the concentration used in the Lipniacki formulation (0.06 M) (19) or nearly one-third of that used in the Hoffmann model (0.1 M) (18). Parameter Sensitivity Analysis-Because the dynamic imaging experiments generate only 15 data points per cell, inferring all 30 parameters in our ODE model from the experimental data for each cell would lead to an ill conditioned inverse problem. Instead, only a subset of the model parameters needs to be inferred. To determine which parameters to infer, a sensitivity analysis was performed to assess the relative influence of parameters affecting the monophasic nuclear NF-B response. To characterize this response quantitatively, we chose four characteristic output quantities: the stationary value for the unstimulated period (Z 1 ), the timing of the maximum value (Z 2 ), the maximum value itself (Z 3 ), and the final value (Z 4 ) (Fig.  5A). For the sensitivity analysis, these 30 parameters were perturbed over a range of values that includes the nominal values (19), shown in Fig. 5B (see supplemental material for a discussion of the sensitivity analysis method used). Besides the three parameters describing the initial state and geometry (k v , T c , and initial concentration of NF-B), the rate of degradation of the bound form of the IB␣ protein in the IB␣⅐NF-B complex (c 6a ), the rate of IB␣ protein production (c 4a ), the IB␣ mRNA degradation rate (c 3a ), and the rate of IB␣-inducible mRNA synthesis (c 1a ) showed the strongest influence on the dynamic behavior of the system. Given the strongly correlated effect of c 4a , c 3a , and c 1a , we picked c 4a as a parameter of inference along with, c 6a , k v , T c , and the initial concentration of NF-B⅐IB␣. Furthermore, because the NF-B nuclear import rate (i 1 ) has the next highest sensitivity measure and has the largest influence on the peak timing of the nuclear/total NF-B ratio (Fig.  5B), it was added to the set of six parameters for inference.

TABLE 1 Ratio of cytoplasmic to nuclear volume (k v ), measured from sequential axial images from control (C) or stimulated (TNF␣) A549 EGFP-RelA cells
Single Cell Bayesian Inference of Model Parameters for 20 Individual Cells-Parameter inference was performed with the experimental noise characteristics s p and s q for all 20 measured NF-B translocation profiles individually. In the forward model simulations, the initial NF-B⅐IB␣ concentration was set to its corresponding parameter value, and the amount of NF-B in the inactive reservoir was computed to satisfy the parameter T c . Initial concentrations for other NF-B-containing molecular species were set to 0, and concentrations for other molecular species were set to an appropriate initial guess. The simulation was then run for 20 h in an unstimulated regime in order to obtain a steady state that mimics the experimental configuration prior to stimulation. After this initial equilibration, stimulation with TNF␣ was switched on in the model, and the simulation was continued for the duration of the experiments. Fig.  6A shows the experimental data and the ODE model prediction  A, all quantities are related to the nuclear/total NF-B ratio values; Z 1 is the stationary value of the unstimulated period, Z 2 is the time of the peak, Z 3 is the peak value itself, and Z 4 is the final value. B, an overall sensitivity measure S i ϭ ͉␣ 1i ͉ ϩ ͉␣ 2i ͉ ϩ ͉␣ 3i ͉ ϩ ͉␣ 4i ͉ is plotted, where ␣ ki is the sensitivity of the kth output quantity of interest with respect to the ith input parameter, where k ϭ 1, 2, 3, and 4, and i ϭ 1, 2, . . . 30. OCTOBER 28, 2011 • VOLUME 286 • NUMBER 43

JOURNAL OF BIOLOGICAL CHEMISTRY 37749
with the maximum a posteriori (MAP) inferred parameter values for one sample cell, cell 8. The figure also shows the model prediction based on the mean values of the parameters, together with the associated posterior predictive uncertainty that incorporates the experimental error with the uncertainty associated with the finite number of data (36). The corresponding MCMC chains for six inferred model parameters are shown (Fig. 6B).
To quantify the uncertainty in and the correlations between the inferred parameters, we analyzed the joint posterior distributions of the six inferred parameters for the given cells (Fig.  6C). All of the conclusions below hold in general for all cells, although they are only illustrated on the single, eighth cell. The posterior distribution summarizes the current state of knowledge about all of the inferred parameters based on the available data. These posteriors were estimated from the last 40,000 samples of the corresponding MCMC chain (Fig. 6B illustrates the full chain) to ensure that the samples were drawn from a chain that had reached a stationary state.
The posterior distributions incorporate the information from the prior distributions and the likelihood function. In the absence of any quantitative information on rate constants, the prior distributions typically provide only reasonable physical bounds on the values of these parameters. The likelihood function, on the other hand, is based on fitting the model to the experimental data; hence, generally the posterior values are driven by the likelihood function. This explains, for example, the strong positive posterior correlations between the total NF-B concentration T c and initial IB␣⅐NF-B concentration (Init) (Fig. 6C). Indeed, as the total concentration increases, in order for the model to fit the data comparably well, the model should use higher initial concentration because the data are based on ratios of concentrations only. Similarly, we notice high positive posterior correlations between c 4a and c 6a , reflecting the fact that IB␣ mRNA translation and IB␣ degradation from the IB␣⅐NF-B complex have opposite effects on the nuclear/total NF-B ratio (Fig. 6C). Increasing c 6a tends to increase the amount of translocated NF-B by accelerating the translocation, whereas increasing c 4a reduces the amount of translocation of NF-B by accelerating the negative feedback through IB␣, discussed below.
The single cell analysis produced good fits with the experimental data based on 80,000 MCMC steps for each of the 20 cells (Fig. 6D). The shaded region corresponding to the posterior predictive uncertainty on the inferred model predictions for each of the 20 cells suggests that the model has been able to capture most of the inherent variation of the experimental data. The MAP estimates for the inferred parameters are shown in Table 2 and Fig. 6E together with the corresponding mean estimates and the associated error bars. The greatest variation in the MAP values for the six inferred parameters across these 20 cells is seen for the c 4a and c 6a rate constants. Note that the inferred MAP values for k v and T c are very similar for all cells and are close to the maxima of their corresponding prior distributions. This indicates that the inference is driven mainly by the initial information (i.e. the priors for these parameters). Namely, the MCMC chain tends to move the MAP estimate toward the maximum of the prior distribution if the corresponding parameter does not have much effect on the likelihood function. Obviously, although the ODE model is able to  reproduce the translocation dynamics observed in each of the 20 cells, some of the inferred parameters may not be biologically realistic. One possible reason for obtaining unrealistic parameters in the inference is the limited amount of experimental data available to constrain the inferred parameters. The ratio of experimental data versus the number of parameters to be inferred can be improved by recognizing that some parameters, such as reaction rate constants, can be expected to be the same in all cells. We therefore performed an inference of the ODE model parameters on all 20 cells simultaneously. Multicell Bayesian Inference of Model Parameters Using Data from 20 Cells Simultaneously-For the inference of parameters in a single ODE model based on multiple-cell NF-B translocation profiles simultaneously, the reaction rate constants (c 4a , c 6a , and i 1 ) were assumed to be the same among all 20 cells, whereas each cell was allowed its own, cell-specific values for the volume ratio k v , the initial NF-B⅐IB␣ concentration, and the total NF-B concentration T c . We believe these assumptions are justified because the biochemical processes governing the IB␣ translation (c 4a ), the IB␣ proteolysis (c 6a ), and NF-B nuclear import rate (i 1 ), are unlikely to change from cell to cell, whereas our experimental data indicate that the NF-B concentrations and cell volumes do vary across the cell population. Thus, with a total of 20 cells, there were 63 (3 ϫ 20 ϩ 3) parameters to infer. Note that although this presents a higher dimensional inference problem than in the previous case of single cell model inferences, it is better constrained because in the single cell inferences, a total of 120 (6 ϫ 20) parameters had to be inferred against the same amount of experimental data. As shown in Fig. 7A, predictions with the inferred parameters for this multiple cell model show good agreement with the respective experimental data for 16 of the 20 cells. Also, it is for these 16 cells only that the multicell inference leads to fits that are similar to the inferred single cell model predictions for the corresponding cells (Fig. 6D). Predictions with the globally inferred model parameters for cells 8, 13, 14, and 20 deviated from the experimental data considerably. For these cells, experimental data suggest a much faster rate of nuclear-to-total NF-B decay than indicated by the inferred model. Faster or slower decaying signals compared with the regular kinetics are reflected in the inferred values of c 4a and c 6a . The reaction rate c 4a implies faster appearance of IB␣ in the system due to translation, thereby inactivating the active NF-B signal, which can be seen as quick decay from the translocation peak in four cells (cell 8, 13, 14, and 20) compared with others.
The MAP estimates for the parameters inferred with this multicell analysis are listed in Table 3 and Fig. 7B, together with the mean estimates and the corresponding error bars. Comparing Fig. 6D and Fig. 7A, it is clear that the inferred parameters show a lot more cell-to-cell variability in the multicell inference than in the case where each cell's dynamics were fit individually. Also, smaller error bars for the multicell case are explained by the fact that the multicell inference is better constrained than the single cell inferences (i.e. it has more data available per inferred parameter or, in other words, fewer degrees of freedom to fit the data).
The greater cell-to-cell variation is expected for the inferred parameters because forcing the reaction rates to be shared between all cells now only leaves three parameters instead of six to account for cell-to-cell variability in the observed translocation dynamics. In this context, it is illustrative to study the role of the prior distribution in the Bayesian inference. As shown in Fig. 6E (panel corresponding to k v ), the inferred values for k v in the single cell inference are between 2 and 3 for all cells (i.e. very close to the maximum of the prior distribution for k v ) (supplemental material). This means that in this case, the prior performs a regularizing role in the inference. Because there are not many data to constrain the parameters, the inference problem is borderline underdetermined, meaning that there is not enough information in the data to properly infer all parameters. The Bayesian inference in this case tends to select values for k v that follow the prior distribution (i.e. prior information is used to constrain the parameter values to make up for the lack of information in the observational data). For the multicell inference (Fig. 7B), most of the inferred values for k v appear near the edges of the prior distribution for k v . In this case, there is a higher ratio of data points versus parameters to be inferred, and the posterior distribution is determined more by the likelihood function than by the prior. In order to fit cell-to-cell variations in the observational data, the inferred parameter values also vary from cell to cell. The prior distribution in this case performs more of a bounding role, keeping the inferred parameters from moving outside of the range of feasible values observed in the z-stacking experiments that the prior for k v was derived from. The fact that almost all inferred k v values are near the edges of the k v prior distribution suggests that biologically unrealistic parameter k v values would be needed in order to generate a better fit with the observed data, which is an indication that the inferred ODE model may not be the most suitable for the full set of 20 cells or that other parameters, besides the ones inferred in this analysis, also need to be retuned.
Inference of Model Parameters Considering Two Subpopulations-The above multicell inference analysis suggests that there are two subpopulations in the cells studied here. Notably, four cells (8, 13, 14, and 20) showed fast decay in the NF-B translocation signal, with a clear inflection point in the decay signal, whereas the rest of the cells exhibit relatively slower decay dynamics. To explore this possibility, the multicell Bayesian inference was repeated, but this time for the fast decay subpopulation (cells 8, 13, 14, and 20) and the other, slow decay, subpopulation separately. The results, summarized in Fig. 8A, show much better fits with the data for all cells in the subpopulation-based analysis than in the full 20-cell multicell inference. Their mean estimates are presented in Fig. 8B.
Our experimental data indicated the existence of a reservoir of translocation-competent NF-B that represented only a fraction of the total NF-B. We therefore calculated the fraction of the translocation-competent NF-B given as the abundance of the NF-B⅐IB␣ complex for each of the 20 individual cells based on the parameter estimates. These distributions are plotted in Fig. 9 and agree well with our experimental observations (cf. Figs. 1A and 2 (A and F)).
Simulations of c4a Variations-Our analysis of the multicellular parameters suggests that dynamics in the regulator feedback loop, such as cellular variations in IB␣ translation rate (c 4a ), can affect the decay rate of the monophasic NF-B profile.
To illustrate this prediction, simulations using variations in c 4a were performed (Fig. 10). With increasing IB␣ translation rate, more IB␣ protein is synthesized per unit time than nominally. The newly synthesized IB␣ saturates free NF-B and traps it into the IB␣⅐NF-B complex, reducing free NF-B available for nuclear translocation. As a result, the initial peak of nuclear NF-B is reduced, and despite the presence of a tonic TNF signal, so much IB␣ is produced that the secondary NF-B response can no longer be seen. Conversely, as c 4a is reduced relative to nominal values, more free NF-B is available to translocate, and less IB␣ is present, permitting a secondary translocation response (gray to orange to blue curves). These data suggest that differences between these subpopulations can be attributed to the dynamics of the regulatory feedback loops and are a good objective for further studies.

DISCUSSION
Despite the theoretical treatment of the NF-B signaling pathway, there are few data available for a rigorous inference of the model parameters or to account for major sources of cellular variability in NF-B pathway dynamics. In this study, we applied a modified model of the canonical NF-B pathway that accounts for a distinct subpopulation of translocation-competent NF-B. We applied Bayesian inference to estimate important parameters using time-resolved NF-B translocation measurements at a single cell level. Bayesian inference is well suited for this approach, being able to tolerate large noise levels and to provide a joint probability distribution for correlated parameters. Experimental measurements of cytoplasmic and nuclear volumes and fractional NF-B translocation were used to inform the prior distribution of the cellular geometry parameter k v and the translocation-competent NF-B⅐IB␣ complex. This Bayesian inference allowed a better understanding of the single cell dynamics and insight into the major sources of cellular variability in NF-B signaling.
There are several new insights generated from our study. First, we note that the total abundance of RelA varies considerably from cell to cell, a finding that may contribute to significant cell-to-cell variation in the activity profile of the pathway. This is observed both in our stable EGFP-RelA transfectants and for endogenous RelA in non-transfected cells, indicating external random effects on the NF-B pathway in each cell. This point is important because models produced based on population-averaged experimental data will not accurately capture the variation of individual cellular responses.
Another significant finding from this study is that a large proportion of cytoplasmic NF-B is not induced by a saturating concentration of TNF␣ to translocate into the nucleus. This result is quite surprising because TNF␣ is one of the most potent NF-B activators currently known (25). We provide data that indicate the inactive complex is due to association with the NF-B1 and -2 precursors, cytoplasmic ankryin repeat-containing proteins that sequester a greater molar amount of RelA than that associated with individual IB isoforms. As a result of these experimental measurements, we modified the deterministic model to include two forms of NF-B: one as an inactive reservoir in the cytoplasm and the second as a translocation-competent reservoir as an IB␣⅐NF-B complex. Our inclusion of NF-B reservoir in the model resulted in dramatically better fitting of the model with experimental single cell data. More work will be required to understand how the reservoir exchanges with the translocation-competent IB␣⅐NF-B pool and how stimuli that induce NF-B1/2 processing affect the dynamics of the pathway.
Although it has been generally assumed that the cellular geometry affects the dynamic NF-B response, this critical parameter has been assumed in previous modeling studies. k v determines the concentration of nuclear NF-B as it undergoes nuclear import, with the initial NF-B concentration determining the kinetics of activation of the negative feedback loops. Our experimental measurements indicate that k v in A549 cells is distributed between 2 and 3, values that are lower than previous modeling assumptions. Also, k v does not dramatically change as a result of pathway activation.
Bayesian inference was then applied to this refined NF-B signaling pathway to estimate parameters using the single cell imaging profiles. The choice of the six inferred parameters in the Bayesian methodology was dictated by a global sensitivity analysis of all 30 input parameters around the nominal values. It is important to note that the inferred concentrations of active NF-B (i.e. in the form of the IB␣⅐NF-B complex) are generally about half the concentration used in previous models. Incorporating prior measurement data on k v and the total NF-B concentrations in the inference helped to constrain the inferred parameters to biologically realistic values and made the inference process better conditioned overall. The probabilistic nature of the Bayesian inference allowed us to develop the model on a single cell basis and assess the uncertainty in the inferred parameters.
Using inferred values of k v and the total NF-B concentration (T c ) for all 20 studied cells, we estimated the fraction of reservoir inactive NF-B concentration in these single cells. Fig. 9 clearly shows a high fraction of reservoir NF-B in all cells as well as considerable cell-to-cell variability in this quantity. This may be one of the major sources of cell-to-cell variability because different cells have different amounts of total NF-B and different fractions of initial active NF-B, affecting the dynamics of NF-B translocation and the downstream signaling initiation process.
Although the analysis above suggests several crucial parameters influencing NF-B translocation dynamics that exhibit cell-to-cell variability, those parameters were not sufficient to capture the cell-to-cell variability observed in the present data set. A Bayesian inference of a single ODE model against data from all 20 cells simultaneously, assuming the rate constants are shared between all of the cells, whereas initial and geometric properties of the cells are cell-specific, was unable to fit all cells. This multicell analysis clearly suggested the existence of two subpopulations of cells; four cells showed much faster decay in the nuclear/total NF-B ratio dynamics than the rest of the cells. Separate multicell inference runs for these two subpopulations resulted in much improved model fits with the experi-mental data for all cells. With all other parameters statistically similar for the cells in both subpopulations, both c 4a and c 6a were significantly larger for the cells in the fast decay subpopulation than in the slow decay subpopulation. The observed subpopulations therefore point at another mechanism of cell-to-cell variability. Besides the unlikely pos- sibility that there are two distinct cellular subpopulations with very different NF-B⅐IB␣ degradation and IB␣ translation rates, another hypothesis is that other factors may drive cells into either a slow or a fast decay mode. One possibility is that other parameters coupled to IB␣ degradation or resynthesis rates are under stochastic expression behavior (28), a feature not captured in the current ODE model. For example, stochastic influences may affect IB␣ mRNA generation and consequently change the dynamics of the IB␣ negative feedback loop, leading to cell-to-cell variability in the NF-B translocation dynamics in a pool of otherwise identical cells. Alternatively, other parameters may need to be inferred, or the current ODE model may need to be modified to capture missing dynamics in the IB␣ degradation. More focused experimentation, outside of the scope of the present study, will be required to identify and understand the source of the pathway deviation in the observed subpopulations.
In summary, this study reported and validated novel experimental observations indicating both active and inactive NF-B in cells. The study also performed a novel Bayesian inference of ODE model parameters against single live cell translocation data, using data from experimentally measured cellular compartment volumes and NF-B concentrations as prior information. We observed cell-to-cell differences in inferred values of compartment volumes, initial concentrations of NF-B, and the translation rate of IB␣ along with the rate of degradation of IB␣ from the NF-B⅐IB␣ complex. All of these factors influence the overall NF-B dynamics of cells and are sources of cell-to-cell variability in the observed NF-B translocation dynamics.  The c 4a parameter represents the IB␣ translation rate in the system, which affects the NF-B activity. A simulation of NF-B signaling was performed using a parameter multiplier for the c 4a rate, ranging from 0.000976 to 256 (i.e. 1/512x to 512x, where x is the nominal value of c 4a ). The nominal value used here is 0.5. We plotted the calculated NF-B activity (ratio of nuclear/total NF-B concentration) versus time. Increasing the c 4a rate results in faster decay of the nuclear NF-B activity.