Abstract
Background
Activation of the NFκB transcription factor and its associated gene expression in microglia is a key component in the response to brain injury. Its activation is dynamic and is part of a network of biochemical species with multiple feedback regulatory mechanisms. Mathematical modeling, which has been instrumental for understanding the NFκB response in other cell types, offers a valuable tool to investigate the regulation of NFκB activation in microglia at a systems level.
Results
We quantify the dynamic response of NFκB activation and activation of the upstream kinase IKK using ELISA measurements of a microglial cell line following treatment with the proinflammatory cytokine TNFα. A new mathematical model is developed based on these data sets using a modular procedure that exploits the feedback structure of the network. We show that the new model requires previously unmodeled dynamics involved in the stimulusinduced degradation of the inhibitor IκBα in order to properly describe microglial NFκB activation in a statistically consistent manner. This suggests a more prominent role for the ubiquitinproteasome system in regulating the activation of NFκB to inflammatory stimuli. We also find that the introduction of nonlinearities in the kinetics of IKK activation and inactivation is essential for proper characterization of transient IKK activity and corresponds to known biological mechanisms. Numerical analyses of the model highlight key regulators of the microglial NFκB response, as well as those governing IKK activation. Results illustrate the dynamic regulatory mechanisms and the robust yet fragile nature of the negative feedback regulated network.
Conclusions
We have developed a new mathematical model that incorporates previously unmodeled dynamics to characterize the dynamic response of the NFκB signaling network in microglia. This model is the first of its kind for microglia and provides a tool for the quantitative, systems level study the dynamic cellular response to inflammatory stimuli.
Background
The nuclear factorκB (NFκB) transcription factor is ubiquitously expressed in mamallian cells and regulates the expression of many target genes. In the nervous system NFκB is known to play a key role in the immune and injury responses and in governing normal brain function [1]. During cerebral ischemia NFκB is a primary regulator of the inflammatory response to ischemic injury, affecting cell death and survival [2]. Microglia, the resident immune cells in the brain, are activated following ischemia and play a controversial role in this decision. Microglia respond to injury in part by releasing both cytoprotective and cytotoxic signaling molecules to surrounding cells, many of which are regulated by NFκB [3]. As the dynamics of NFκB activation control gene expression [46], characterizing the dynamics of NFκB activation in microglia is of great interest.
Members of the NFκB family of transcription factors are found in their inactive state as dimers bound to their IkB inhibitor proteins. Upon stimulation by a diverse set of stimuli, NFκB is freed from its inhibitor to coordinate gene expression in a highly specific and tightly regulated manner. The IκBα inhibitor and p65(RelA):p50 NFκB heterodimer are the most extensively studied members of their respective families, and their response to extracellular stimuli illustrates the canonical pathway of NFκB activation (Figure 1).
Figure 1. The canonical NFκB activation pathway. Binding of TNFα trimers to TNFR receptors initiates the canonical signaling pathway by activating the upstream kinase IKK. IKK phosphorylates the IκB inhibitor that is bound to NFκB in the resting state. This targets IκB proteins for the ubiquitinproteasome system, which leads to IκB destruction by the 26S proteasome and release of NFκB. Free NFκB enters the nucleus and activates gene expression of many target genes and induces negative feedback regulation by synthesizing IκB and A20. IκB proteins inhibit NFκB activity by sequestering NFκB from the nucleus to form an inner feedback loop, while A20 attenuates stimulus induced IKK kinase activation further upstream in an outer negative feedback loop.
In the canonical pathway, binding of extracellular TNFα trimers to TNFR1 receptors at the cell membrane initiates NFκB activation. The ligandreceptor complex interacts with several adapter proteins, including TNF receptorassociated factor 2 (TRAF2) and receptorinteracting protein1 (RIP1), which are essential for recruitment and activation of the IκB kinase complex (IKK) [7]. The IKK complex involved in canonical NFκB activation is composed primarily of the regulatory subunit IKKγ (NEMO) and two catalytic subunits: IKKα/IKK1 and IKKβ/IKK2. Upstream signals activate IKK by phosphorylation of the kinase domain of IKKβ, which in turn phosphorylates IκBα on serines 32 and 36 [8]. Phosphorylated IκBα is recognized by the βTrCP containing Skp1CulinRoc1/RBx1/Hrt1Fbox (SCF) E3 ubiquitin ligase complex (SCFβTrCP), which facilitates K48linked polyubiquitination of IκBα and targets it for degradation by the 26S proteasome [9,10].
NFκB is released following proteasomal degradation of IκBα [11] and translocates to the nucleus, where it activates gene expression. Of the hundreds of genes targeted by NFκB [12], two in particular are ikba and a20. The expression of these genes is rapidly induced by NFκB and triggers the synthesis of de novo IκBα and A20 proteins. Newly synthesized IκBα sequesters NFκB from the nucleus to inhibit further transcriptional activity, forming a strong negative feedback regulatory mechanism. The synthesis of A20 proteins creates a second negative feedback loop by regulating the ubiquitination of adapter proteins responsible for activating the IKK complex, thus inhibiting further NFκB activation [13]. Many characteristics that define TNFα induced NFκB activation also underlie cellular responses to many other stimuli, necessitating a thorough understanding of this pathway.
Given the dynamic nature of NFκB signaling and its regulation involving multiple feedback loops, it is necessary to consider the network as a whole when studying this system. The seminal work by Hoffmann and colleagues [4], in which simulation predictions were used in coordination with experimental studies of IκB knockout cells to reveal functional differences among three IκB isoforms, established mathematical modeling as a vital tool for studying NFκB signaling at a systems level. Subsequently a number of researchers have used modeling to investigate various aspects of NFκB activity [5,6],[1418].
Here we develop a mathematical model to describe NFκB signaling in microglia. Beginning with a recently published model structure shown to be capable of predicting NFκB signaling in other cell types [14], we attempt to identify model parameters to match experimental data sets of NFκB and IKK activation obtained from a microglial cell line. The inability of the original model to recapitulate NFκB activation that is consistent with experimental data  regardless of model parameter choice  leads us to expand the model to incorporate previously unmodeled dynamics of the IκBα ubiquitinproteasome degradation pathway. We also find that IKK activation in microglia is highly nonlinear, which prompts refinement of the upstream signaling module. We use the new model to predict the levels of another network component, total IκBα, and are able to validate this prediction experimentally. The results offer a validated model that can be used as a new tool to study the dynamics of NFκB activation in microglia. While we find that many key features of canonical NFκB activation are shared in microglia, the model suggests a potentially more prominent role for the ubiquitin system in regulating the dynamics of NFκB activation. We use numerical analyses of this model to gain insight into how microglia regulate both IKK and NFκB activity in response to inflammatory stimuli. Our sensitivity anlayses emphasizes the dynamic nature of how key system responses are regulated, a feature that may not be apparent from similar analyses. The analysis further highlights the robust yet fragile nature of the NFκB signaling pathway due to the multiple layers of feedback regulation.
Results
TNFα stimulates dynamic NFκB and IKK activation in BV2 microglia
To characterize the dynamics of canonical NFκB activation in microglia, cells from the microglial cell line BV2 were cultured and treated with 10 ng/ml TNFα. Whole cell extracts were collected in triplicate over a time course following stimulation in five identical experiments conducted on different days. ELISA measurements of NFκB p65 DNA binding activity show that NFκB activation in BV2 microglia is strongly induced by TNFα (Figure 2A). Five minutes following TNFα treatment NFκB activation remains near basal levels but increases rapidly thereafter, reaching maximal activity near 20 min. Following the initial peak, NFκB activity declines until approximately 90 min when it returns to a second, smaller amplitude peak. The biphasic NFκB activity profile in BV2 microglia is consistent with the NFκB response to sustained TNFα stimulation observed from population level measurements in many other cell types [19].
Figure 2. Dynamics of NFκB and IKK activation in BV2 microglia treated with TNFα. ELISA measurements of (A) NFκB p65 DNA binding activity and (B) IKKβ kinase activity following continuous stimulation by 10 ng/ml TNFα. Data markers at each time point are sample averages from independent experiments performed on separate days. Error bars indicate one standard deviation of the samples.
To better characterize the inflammatory response in microglia we additionally examined the activation of the upstream IκB kinase (IKK) experimentally. The time course of IKK activity was measured for the first 30 min following 10 ng/ml TNFα treatment in three identical experiments. IKK is rapidly activated, reaching peak levels near 5 min. By 10 min IKK activity sharply drops to below halfmaximal levels and gradually declines to near basal levels over the next 20 min (Figure 2B). This transient profile resembles IKK activation characteristic of the response in most other cell types to high TNFα doses, in which IKK activity peaks between 515 min and drops below 25% of its maximal value by 30 min [20,21]. However, the rapid decline from maximum activity at 5 min to ~33% activity by 10 min is particularly prominent in microglia.
Intermediate steps in the IKKinduced IκBα degradation pathway reconcile the mathematical model with NFκB activation in microglia
Next we sought to quantitatively describe microglial NFκB activation using a mathematical model. While a number of mathematical models for NFκB have been published in recent years (reviewed in [19]), our preference was to begin with a simple description that still captures the essential components of the network. For this purpose we selected a deterministic, ordinary differential equation (ODE) model structure recently published by Ashall et al [14], which was based primarily on an earlier model by Lipniacki et al [18]. This model includes the core architecture of the canonical signaling pathway and was able to predict many key features of NFκB activation in different cell types under a variety of conditions.
We first attempted to identify parameters for the existing model structure to fit the experimental NFκB and IKK activation profiles of microglia. An optimizationbased parameter estimation algorithm was run using many randomly selected parameter values from the parameter space as initial guesses. However, no parameter sets were found that matched microglial IKK and NFκB activity. In particular, the model was unable to qualitatively reproduce the rapid induction and attenuation of IKK activity observed in microglia for any of the parameter sets tested, and NFκB activation was predicted to occur more rapidly than the 5 min delay observed in Figure 2A. The discrepancies between the model and data prompted us to investigate the time interval immediately following TNFα stimulus.
Sensitivity analyses were performed on the model to quantify the relative contributions of each of the system parameters to the concentration of free NFκB during the first 10 min given the large mismatches between the model and data in this interval. Only seven of the original 26 system parameters have appreciable effects on NFκB activity during this time based on their timeaveraged sensitivity scores (Figure 3A, Additional file 1: Figure S1). Notably, the most significant parameters comprise the rates governing IKK activity, IKKinduced phosphorylation and degradation of bound IκBα, nuclear import of NFκB and its association with IκBα, and the ratio between the volumes of the cytoplasm and nucleus. No parameters governing transcriptional regulation or other downstream processes have significant effects on NFκB activation during this early time interval as evidenced by their very small sensitivity scores. Moreover, this ruled out the possibility that feedback from other IκB isoforms (e.g. IκBε) not included in this model could be added to account for the discrepancies in the dynamics. This suggested that the brief delay in the initiation of NFκB activation observed in microglia was likely due to unmodeled dynamics involved in the IKKdependent degradation of IκBα or to dynamics in the upstream signaling pathway governing IKK activation, allowing us to restrict our initial attention to only a subset of key upstream parameters.
Figure 3. New model structure required to characterize NFκB activation in microglia. (A) NFκB activity during the first 10 minutes following stimulation was only highly sensitive to seven of the 26 rate parameters. (B) By using an IKK signal derived from experimental measurements as the model input, the outer feedback loop can be removed (indicated by gray lines), isolating the downstream NFκB activation module with IκBα feedback. Similarly, once the concentration of nuclear NFκB is known, this signal can be used to drive the upstream IKK activation network independently of the downstream module. (C) Model structure from the original model (top) and the new model (bottom). (D) Simulations with parameters estimated for the existing model (dashed line) and the new model (solid line) using the experimental IKK curve as input. The inset provides a detailed view of the model fits during the initial activation phase. (E) The results of 1980 randomly initialized parameter estimates for each model were checked for statistical consistency with the data using Fisher's Method (see Methods) and binned according to pvalue. No estimated parameter sets with the original model achieved a Pvalue >10^{7 }(red), while nearly half the estimated parameter sets with the new model (blue) had P > 0.01.
Additional file 1. Supplementary text and figures. The pdf contains supplementary text describing development of the mathematical model; Tables S1S3 which list the model species, reactions, and rate parameters; and Figures S1S8 that provide more detailed simulation results.
Format: PDF Size: 702KB Download file
This file can be viewed with: Adobe Acrobat Reader
To more easily explore these possibilities and to facilitate model development, we first considered the downstream network independently of the upstream IKK activation network. IKK interacts with the downstream module only through its enzymatic phosphorylation of IκBα and through feedback inhibition from A20 (Figure 3B). We isolated the downstream network by breaking the outer A20 feedback loop and using the interpolated experimental IKK activation data as the model input in a manner resembling previous work by others [22].
With the IKK profile fixed as the model input, the least squares parameter estimation procedure was repeated with certain parameter values and biological features constrained by the literature (Additional file 1: Table S2). Simulations of the existing downstream model with the estimated parameters predicted free NFκB levels increasing sooner than what was detected in microglia, as was also the case for the full model (Figure 3D). To test whether this result was limited to a particular set of values or held more generally, many additional estimates were obtained starting from initial values randomly sampled from the parameter space using both a leastsquares objective function and an alternative objective function adapted from the parameter estimation method proposed by [23]. Following the methodology in [23], we applied an a posteriori statistical test based on Fisher's Method to check whether model simulations at each estimated parameter set were consistent with the experimental data, taking into account measurement errors in the data (see Methods). The results showed that with the original model structure, 100% of the estimated parameter sets had Pvalues < 10^{7 }(Figure 3E), leading us to conclude that the original model could not produce dynamics consistent with the data.
Taken together with the sensitivity results showing that very few system parameters significantly affect NFκB activation during the first 10 min of activation, this strongly suggested there were likely unmodeled dynamics within the IKKinduced IκBα degradation pathway. We next investigated whether the model could be modified in a biologically meaningful way to incorporate missing dynamics and to better fit the data.
The original model structure describes IKKdependent IκBα degradation in two steps: phosphorylation of IκBα catalyzed by IKK, and degradation of phosphorylated IκBα (Figure 3C). However, this twostep description omits many intermediate steps which occur prior to IκB degradation by the 26S proteasome [7]. We therefore extended the model to include two intermediate reactions following IκBα phosphorylation and preceding IκBα degradation, which we posited might be sufficient to account for the missing dynamics. The reactions roughly correspond to recognition of phosphorylated IκBα by an E3 ligase intermediate, and attachment of a ubiquitin (Ub) chain to the substrate (Figure 3C). It must be noted that each of these reactions potentially encompasses numerous intermediate steps and may not correspond directly to the reactions as they are described here; however, the mechanistic details of this pathway obtained from the literature provide a biological basis for developing this model.
With the new model structure in place, the parameters corresponding to the new stimulusinduced IκBα degradation reactions were estimated using the optimization algorithm while fixing all other parameters downstream of IκBα degradation to their previously estimated values. Remarkably, parameters were found to closely match microglial NFκB activation, decreasing the data fitting error by nearly 67%, with over a 9fold improvement during the first 20 min in particular (Additional file 1: Figure S2). Reestimating the other parameters with the modified model provided even better agreement with the data, further reducing the fitting error from 0.67 to 0.30 (Figure 3D). The consistency between simulations of the new model and the data was assessed using the a posteriori statistical test as before. At these parameters the test yielded a Pvalue of 0.038, implying that the null hypothesis could not be rejected with a high significance level. This result was corroborated by obtaining a large number of parameter estimates and finding that nearly 50% of the estimates with this model structure had P > 0.01 (Figure 3E).
These results provide strong evidence that the addition of dynamics roughly corresponding to the steps involving phosphorylated IκBα recognition and binding by the E3 ligase, polyubiquitination, and proteasomal degradation is sufficient to account for the slightly delayed NFκB activation observed in microglia.
Nonlinearities in IKK activation and inactivation produce the rapid transient IKK activity in microglia
We next focused our attention on the upstream signaling pathway governing IKK activation in response to TNFα stimulation. The upstream signaling module was decoupled from the downstream model by using the concentration of free nuclear NFκB produced by the downstream module as a fixed model input (Figure 3B). This enabled us to consider only the reactions immediately governing IKK activity and its regulation by A20, again greatly simplifying the model development task.
The original upstream model in [14], which includes IKK cycling among three states (native, active, and inactive) and feedback from A20, was unable to adequately fit either the rapid activation or deactivation of microglial activation (Figure 4B). Therefore, we examined ways in which the model could be modified consistent with the biology to better correspond with the data.
Figure 4. Development of the upstream model and full model fits. (A) The upstream model with A20 inhibition at two points was modified to include nonlinear IKK activation and inactivation rates, indicated by dashed lines. (B) Simulations with the newly developed upstream model (solid) provided excellent agreement with the model (P = 0.854, SSE = 0.038), much better than what was possible with the existing model structure (dashed) (P = 0.0002, SSE = 0.661). With the newly developed upstream and downstream modules integrated into the full model, new parameter estimates were able to fit both NFκB (C) and IKK (D) activity in microglia. (E) Model predictions for the total amount of IκBα (solid line) are compared with experimental measurements of total IκBα protein in BV2cells following 10 ng/ml TNFα treatment.
Activation of the IKK complex at the biomolecular level involves the recruitment and assembly of a signaling complex following TNFα binding to its receptor, as well as numerous posttranslational modifications to the complex subunits before IKK is activated by phosphorylation at two residues within its kinase domain [7]. Although other studies have attempted to model the upstream pathway in a greater level of detail [2426], many of the details are still being resolved and we opted to retain the basic IKK cycling description from [14]. The activation reaction rate was changed from a linear function to a nonlinear Hill equation as a coarse approximation to the many intermediate steps involved in IKK activation.
The quick attenuation of IKK activity following its induction is essential to proper signaling and the resulting biphasic NFκB activity [20,27]. IKK reportedly undergoes hyperphosphorylation at 9 or 10 residues in the Cterminal, which was found to significantly decrease kinase activity in cells [27]. We posited that potential cooperativity in IKK inactivation due to autophosphorylation may lead to nonlinearites in the inactivation rate equation of the model. Accordingly the linear reaction rate was changed to a nonlinear Hill equation.
Feedback from A20 in the published model was proposed to inhibit the transition of inactivated IKK back to its native state [14]. Because we were unaware of any biological basis for such a mechanism, we adopted two mechanisms of A20 interaction (Figure 4A) that had been identified in the literature and had also been included in prior models. The first is direct inactivation of the IKK complex by A20 protein, a mechanism reported in [28,29] and previously modeled in [18]. We used the identical mathematical description of this interaction from [18] in our model. Secondly A20 is known to inhibit activation indirectly through its ubiquitinediting activities of upstream signaling components [13]. This mechanism has been included in previous models that have a more detailed description of the upstream signaling pathway [25,26]. We adapted this second interaction to our model by assuming that A20 attenuates the rate of TNFinduced IKK activation in a concentration dependent manner.
Parameter estimation was performed using the newly developed upstream model with fixed nuclear NFκB as the model input. Parameters were found for which the model produced excellent agreement with microglial IKK activation (Figure 4B), decreasing the fitting error by more than an order of magnitude compared to the best fit achieved with the original upstream model (error 0.04 compared with 0.66). Consistency of the predictions using the new upstream model with the IKK data is more statistically significant (P = 0.85 compared to P = 0.0002) than with original model structure from [14]. Remarkably, we observed that the best fits with the new model were achieved with high Hill coefficients (>3) for IKK inactivation, suggestive of a highly cooperative mechanism in the underlying biological process (Additional file 1: Figure S3).
The newly developed upstream and downstream signaling modules were integrated to form the full model characterizing both IKK and NFκB activity in response to persistent TNFα stimulus (Additional file 1: Tables S1S3). Model predictions using the parameter sets estimated from the isolated signaling modules, while giving good agreement during the first 30 min, predicted a higher amplitude second phase of NFκB activity (Additional file 1: Figure S4), which was inconsistent with the data (P < 10^{6}). Numerical investigation showed this more oscillatory behavior predicted by the integrated model was due to small changes in the later activation profile of IKK predicted by the upstream model, which had been assumed to remain at a constant, low level when developing the isolated downstream signaling module. After increasing the rate of IκBα nuclear import and reestimating the A20 feedback and IKK recycling rates, the newly developed model was able to provide good agreement with the data, with fitting errors of only 0.34 for NFκB (P = 0.013) and 0.43 for IKK (P = 0.291) (Figure 4C and 4D).
Model prediction validated experimentally
Given that the model was developed using a limited set of data from IKK and NFκB activation, we next sought to test its ability to predict the dynamics of other model species for which no information was used during parameter estimation. The model was first simulated to obtain the levels of total cellular IκBα protein following TNFα stimulus (Figure 4E). The model predicted that the level of protein stays relatively unchanged during the initial delay, but begins a decline by 5 min. At 20 min, the model predicts that IκBα protein levels have been reduced beyond half of their initial amounts.
To test this prediction experimentally, BV2 cells were again treated with 10 ng/ml TNFα, and levels of total cellular IκBα were measured at several time points after treatment using ELISA. The results of the experiments were normalized with respect to the initial quantities and compared with the simulation predictions (Figure 4E). The experimental data were in excellent agreement with the predicted IκBα levels, providing a level of experimental validation to the model.
Model analysis highlights robustness properties of the network and a dynamic role of feedback regulation in both NFκB and IKK signaling
The model was next analyzed using sensitivity analysis to gain deeper insight into how the different components of the system interact to regulate the dynamic NFκB response in microglia. Sensitivity analyses of the NFκB regulatory network have been performed previously [3033], and have provided significant contributions to understanding how the system operates. Here we expand upon these studies by considering the dynamic trajectories of the sensitivity coefficients, and examining how the sensitivity of the system response with respect to network parameters changes with time.
The normalized sensitivity coefficients for NFκB activation were solved and plotted as heat maps to illustrate the dynamic relationship between the signaling components and the system response (Figure 5A).
Figure 5. Model analysis of transient IKK and NFκB activation in microglia. Sensitivity analysis of the model shows the dynamic NFκB response (A) and IKK response (B) are regulated differently by different groups of parameters depending on the time interval of interest. See Additional file 1: Tables S2S3 for complete descriptions of the parameters. (C) Parameter scans were used to find the Euclidean distance between the nominal and perturbed NFκB responses as the parameters varied over four orders of magnitude. Dark colors indicate little change, while lighter colors indicate large changes. (D) Simulated trajectories from the parameter scans for selected parameters from each group show distinctions in how the system is affected by up to +/ 10fold changes from the nominal values. See also Additional file 1: Figures S7 and S8.
The sensitivity results clearly show that the NFκB response is nearly completely insensitive to variations in some rate parameters (rows of light green), but also moderately or highly sensitive to others (dark red and blue) (Figure 5A), consistent with earlier results which found that only a relatively small number of network parameters signifcantly influenced NFκB activity [30]. A notable feature of our analysis is that, with the exception of the NFκB nuclear shuttling rates (ki1 and ke1) for which the sensitivity scores remain high throughout the entire response, NFκB activity exhibits highly dynamic sensitivity with respect to most other parameters. In other words, there is a strong temporal component to the regulation of NFκB activity, where variations in different parameters can exhibit great influence over certain phases of activity but have only marginal effects on activation during other time intervals (Figure 5A and Additional file 1: Figure S6). The first 20 min of NFκB activity is predominantly influenced by the rates for IKKinduced phosphorylation, ubiquitination and degradation and also IKK activation, with little contribution from the feedback parameters. As IκBα is degraded and free NFκB ascends towards its maximal activity, the nuclear shuttling rates of free NFκB have the greatest effect.
However, the system shows extreme sensitivity to rates controlling the inner and outer feedback loops. The system is very senstive to the rates for induced IκBα synthesis and its association with NFκB during a time period coinciding with the decline of the first peak, with synthesis and binding rates negatively affecting NFκB activation. The rate of conversion of inactivated IKK back to native IKK (kp) also is among the most significant parameters in the attenuation of NFκB activity. While NFκB activity is at its lowest levels between 6090 min, the stability of the remaining IκBα transcripts and the induced phosphorylation, ubiquitination and degradation of IκBα exert more influence on free NFκB levels. The second peak of NFκB activity is regulated greatly by the nuclear import rate of free IκBα, as evidenced by the high sensitivity of ki3a only during this time period. Feedback from IκBα again has highly significant contributions to the dynamics of the second peak, with induced synthesis of IκBα and its affinity to unbound IκBα having very high sensitivities.
The NFκB response is also highly sensitive to the outer A20 feedback loop in a timedependent manner. The rates for IKK inactivation by A20 significantly affect the termination of initial NFκB activity as well as the second phase of activity. This effect is actuated through inhibition of the activation of IKK that has recently been converted from the inactive form and made available for activation; feedback from A20 inhibition of IKK activation has a less substantial role on the dynamics in the model. The outer feedback parameters governing A20 act in opposition to the IKK recycling rate (kp) to regulate this response, made clear by the opposite signs of sensitivity values throughout the response.
Although many features of the NFκB response have been studied previously using sensitivity analysis, little attention has been paid to the dynamic sensitivities of IKK. We therefore assessed parameter sensitivities of IKK activation in the same way as just described for NFκB (Figure 5B). IKK activity is sensitive to fewer parameters than NFκB, which is expected due to fewer reactions involved in the upstream module, and its only direct interaction with the downstream signaling pathway occurring through feedback from A20. As with NFκB, the IKK sensitivities are also highly dynamic, emphasizing the dynamic nature of its regulation during the initial transient and late, lowactivity phase. The initial peak only exhibits sensitivity to the activation rate (ka) and inactivation rate parameters controlling the magnitude (ki) and the dissociation constant (kmmi). Twenty minutes after the initial stimulus when IKK is mostly in its inactivated form, the response becomes highly sensitive to the IKK recycling rate and to A20 synthesis, degradation, and negative feedback rates which constitute the outer feedback loop. The late phase IKK response is also relatively sensitive to the rates governing IκBα induced synthesis and transcript stability, and to a lesser extent to its induced degradation of IκBα protein, which indicates that the dynamics of IKK are still highly coupled to the inner feedback loop of IκBα despite the absence of direct crosstalk reactions.
While sensitivity analysis with respect to small variations is informative, the nonlinear nature of the system makes it possible that the results may be different when large magnitude changes to the parameters are considered [30]. Robustness of the system response to large changes in parameter values was therefore assessed by varying each parameter over four orders of magnitude and computing the Euclidean distance between the nominal NFκB response and the NFκB response simulated at these perturbed parameters (Figure 5C and 5D). NFκB activity remains relatively unchanged when many of the parameters for nuclear shuttling and IκBα protein degradation are changed to values which differ substantially from their estimated values, indicating that the system response is relatively robust to changes in these parameters (Figure 5C). Examination of the trajectories at parameter values spanning two orders of magnitude shows that indeed the response remains similar when the protein degradation rates are varied by large amounts, and that altering the nuclear import rate of IκBα changes the amplitude of the second peak but retains an otherwise similar profile (Additional file 1: Figure S7). Consistent with the sensitivity results (Figure 5A) in which NFκB was insensitive to activation and inactivation rates for IKK, the NFκB response is robust to changes in these parameter values (Figure 5C). Only extremely large changes in the IKK activation rate parameters significantly alter the response, with much higher activation rates leading to a more oscillatory response (Additional file 1: Figure S8). The parameter scans also show that the system tolerates up to 5fold changes in the new IκBα induced ubiquitination and degradation parameters while maintaining a similar NFκB response, but with the timing of the first peak slightly shifted (Figure 5D and Additional file 1: Figure S7). Decreasing the rate further, however, decreased the amplitude of the response significantly. Surprisingly the system is relatively robust to the nuclear import and export rates (ki1 and ke1), a result which is unexpected given the sensitivity analysis results in which these rates were among the most sensitive. Large changes in these parameters alter the level of damping in the second phase of the response, but the initial peak remains nearly identical (Figure 5D and Additional file 1: Figure S7).
While the system response is robust to large changes in many of the parameter values, the system is much more responsive to changes in the reaction rates involved in both the inner IκBα and outer A20 feedback loops. In particular, the NFκB activation profile changes significantly when the rates of induced transcription or translation are changed only a small amount, as indicated by the large distance between the nominal and perturbed trajectories at these values (Figure 5C). Changes in these parameters by 3fold significantly alter how quickly the response is attenuated and change the frequency of the second phase of activity (Figure 5D and Additional file 1: Figure S7). Similarly, the distance remains small for only a relatively narrow range of rates near the nominal values for most A20 feedback parameters, indicating that the system response changes appreciably when these rates deviate substantially from their nominal values (Figure 5C). Large changes in the A20 feedback loop parameters significantly alter both the amplitude and timing of the second peak and how quickly the first peak is attenuated, but leave the early dynamics relatively unchanged (Figure 5D and Additional file 1: Figure S8).
Discussion
Our quantitative experimental studies show that microglia share many general features of canonical NFκB activation observed in many other cell types [19]. Namely, microglial NFκB activity exhibits a biphasic profile with a high amplitude first peak followed by a damped loweramplitude second phase (Figure 2). NFκB activation begins following a brief delay of nearly 5 min and reaches a peak near 2025 min, resembling profiles observed in other studies with immortalized mouse embryo fibroblasts (MEFs) (see, e.g. [4,34]). The second phase of activity appears to be lower amplitude and more heavily damped than that observed in fibroblasts [4,20], although differences in experimental measurement techniques make direct comparison difficult. The observed damping may reflect asynchronous and oscillatory responses at the single cell level [5,35]. IKK activity in microglia also resembles the tightly constrained IKK profile in other cell types that consists of a fast initial peak occurring 510 min followed by rapid downregulation to low levels of activity at later times [20]. One distinction of IKK activation in microglia, however, appears to be the severe attenuation of IKK activity by 10 min following stimulation, only 5 min removed from peak activation levels (Figure 2B).
Despite the general similarities in NFκB and IKK activation between microglia and other cell types, a recently published mathematical model of the signaling network [14] was unable to recapitulate the nuances of the rapid attenuation of IKK activity simultaneously with the brief delay in the onset of NFκB activity in microglia. Noting that the largest discrepancies between the data and model simulations occurred within the first 20 min of activation, we used this information together with insight gained from sensitivity analysis to develop a new model that is able to match both IKK and NFκB activity in this cell type.
The new model was developed in a modular fashion, which was made possible by collecting ELISAbased measurements of IKK in addition to measurements of NFκB activity and by exploiting the multiple feedback structure of the network. First the IKK data set from microglia was used to develop the downstream signaling module independently of the outer feedback loop, then the upstream signaling pathway was modified to fit IKK activation data, and finally the two modules were integrated to form the full model for which the parameter estimates were refined. The novel downstream signaling pathway includes additional reactions preceding stimulusinduced IκBα degradation, which are sufficient to capture the delayed onset of NFκB activity observed in microglia (Figure 3D and Figure 4C). The mathematical representation we use to describe the additional dynamics is rather basic, yet captures effects that are likely significant at the biomolecular level. We attribute the intermediate model reactions to key steps in the ubiquitination pathway that implicitly have been lumped together in prior models.
Ubiquitination of IκBα is typically thought to occur almost instantaneously following its phosphorylation by IKK [36]. Consistent with this view, recent in vitro kinetic studies revealed in exquisite detail that the SCFβTrCP E3 ligase sequentially adds ubiquitin (Ub) molecules to phosphorylated substrate to form a polyubiquitin chain able to be recognized by the proteasome in a process lasting only seconds after the first Ub molecule has been added [37]. However, the same study [37] also demonstrated that the addition of the first Ub to the substrate is the rate limiting step and occurs with low efficiency during a single encounter between enzyme and substrate, suggesting that any cellular differences affecting how efficiently the initial Ub is conjugated will contribute appreciably to the dynamics. One such possibility for the differential ubiquitination dynamics is celltype specific expression of the E3 ligase components, such as the Fbox protein, βTrCP, which recognizes phosphorylated IκBα [10]. A smaller pool of βTrCP available to bind IκBα, either as a consequence of reduced expression or increased competition with other substrates such as βcatenin, potentially alters how efficiently the substrate is recognized and hence affects the dynamics. Alternatively differential expression of the two βTrCP isoforms  βTrCP1 and βTrCP2  may in part account for the altered response in microglia, as studies using genetic knockouts of βTrCP1 found that TNFα induced IκBα degradation was impaired but not prohibited [38]. Others have posited that the unstable βTrCP2 isoform may be stabilized by increased levels of phosphorylated substrate [39], allowing the possibility that microglia express βTrCP2 in excess of βTrCP1 and thereby have altered ubiquitination dynamics.
Besides potentially less efficient recognition of IκBα by βTrCP, another possibility is that the normally rapid polyubiquitination of IκBα occurs less efficiently in microglia due to smaller quantities of Nedd8ylated Cul1 in the SCF complex. Conjugation of only a small fraction of Cul1 with Nedd8 greatly increases the efficiency of ubiquitination of IκBα without affecting the association between βTrCP and phosphorylated IκBα [40] due to facilitated recruitment of Ublinked E2 to the E3 complex [41]. It follows then that different levels of Nedd8 or the Nedd8conjugating enzyme, Ubc12, could likely contribute to delayed ubiquitination in microglia. Although we cannot decisively point to a particular mechanism as the source of the additional dynamics needed to match the data in microglia, there are many plausible mechanisms which may warrant further study in the future.
The new model structure indicates a more prominent role of the ubiquitinproteasome system in regulating NFκB activation dynamics, which merits consideration of what are its functional implications on how microglia respond to inflammatory stimuli. Analyses of the model show that the ubiquitinrelated parameters have large effects on the initial activation of NFκB and a relatively smaller role in regulating later dynamics (Figure 5A). Parameter scans validate this, as large changes in these parameters change the timing of the first peak by as much as 15 min and alter the amplitude and timing of the later response somewhat (Figure 5C and 5D, Additional file 1: Figure S7). This suggests that altered ubiquitination signaling may be important to regulating the timing of the initial response, but how this affects gene expression and cellular function is not clear at present.
Substantial modifications to the upstream signaling pathway are required to fit the new model to the microglial IKK activation data. The TNFαinduced IKK activation and inactivation reaction kinetics are changed from first order linear massaction rates to nonlinear Hill equations in the new model. We note that the new model differs from [14] in that it includes mechanisms of A20 feedback that more closely reflect the known biology [13,42], but these mechanisms have also been modeled in previous studies [18,25,26]. The nonlinear reaction rates are essentially black box descriptions of a complex upstream signaling network but allow the model to fit the microglial IKK data remarkably well (Figure 4D). Interestingly, the best agreement with the data is obtained with large Hill coefficients (>3) for the inactivation rate (Additional file 1: Figure S3). This may correspond to cooperativity involved in autophosphorylation at 9 or 10 serines in IKK [27]. Additionally, while autophosphorylation decreases phosphorylation in some cells [27], this effect is not observed in all cells [43], which leaves open the possibility that mechanisms besides autophosphorylation are responsible for the rapid nonlinear deactivation in microglia. Although nonlinearities in the activation and inactivation rates are necessary to match the IKK data well in microglia, they do not appear to have a significant influence on the resulting NFκB activity, as indicated by our parameter scans (Figure 5 and Additional file 1: Figure S8). Similar findings have been reported elsewhere [20,22,26], and suggest that cells respond robustly to TNFα stimulus by producing an initial peak of NFκB activity via transient activation of IKK, even in an uncertain environment in which the precise IKK levels may deviate quantitatively but qualitatively remain the same.
In contrast to the parameters governing initial transient IKK activity, our model analyses indicate that the signaling components which regulate later phase IKK activation also exert significant control over NFκB activation (Figure 5, Additional file 1: Figure S8). Key among these is feedback inhibition by A20, which is known to modulate late phase NFκB activity through its inhibition of IKK activity [25,42]. Our analysis suggests that direct A20 inactivation of IKK contributes more to later regulation than feedback inhibition of IKK activation, although more detailed models are likely to provide better insight into the complex regulatory role of A20. The analysis also shows that the inner feedback loop of IκBα is significant in later regulation, emphasizing the interconnected nature of the system.
The sensitivity analyses of the new model presented here provide new insights into how this signaling pathway is regulated. In particular, we show by examining the temporal evolution of the sensitivities that there is a strong temporal component to system regulation (Figure 5A and 5B). Previous studies have used sensitivity analysis to identify the key parameters affecting the NFκB response. These results have typically been reported by ordering the parameters based on the sensitivity scores observed for certain features of the response like the timing and amplitude of NFκB[30], the L^{2}norm of the dynamic sensitivities [33], or a combination of several dynamic features [32]. While the insights afforded by such analyses are valuable, they can potentially obscure information about the dynamics that are of practical value for model development and parameter estimation. Consider for instance the development of the present model. A reasonable strategy to determine where to modify the model to account for the NFκB delay might be to start by examining reactions described by the most sensitive parameters as suggested by the literature. However, if these sensitivity rankings are based largely on the effects on later dynamics as opposed to the initial activation, as is the case for all of the feedback parameters, then trying to modify these parameters would be much less effective. Thus, results from our dynamic sensitivity analysis can be of particular importance when trying to identify how to modify a model to correct discrepancies between model simulations and data, as it provides valuable information.
It is important to note that our particular model, which is developed to reproduce population average measurements of IKK and NFκB activity in microglia, is not unique and other models are capable of producing the same dynamics. It may be desirable in different contexts to extend or otherwise modify this model to explore aspects not considered here. For instance, delayed negative feedback from the IκBε isoform may also contribute substantially to later phase NFκB signaling dynamics [35,44], but is omitted from the present model. It may be useful to extend the model to include interactions from IκBε in future studies. Using data from bulk population level averages also masks asynchronous NFκB oscillations at the single cell level [5,14,34,45]. Thus a different approach, such as simulating the deterministic model with random parameter distributions [34] or using stochasticdeterministic hybrid models [6,14,26], may be more appropriate when specifically considering individual cell responses.
The analysis from this model for microglial NFκB activation clearly portrays the canonical NFκB response on one hand as very robust: cells are able to parse extracellular signals into transient IKK activation to produce a quick and dynamic rise in NFκB activity, even in the face of uncertainty in many of the reaction rates in both the upstream and downstream pathways. This finding is consistent with sensitivity analysis of related models, in which the response was found to be largely insensitive to the majority of the rate parameters [30]. On the other hand, this analysis reveals the highly responsive nature of the network, evident from the high sensitivity and low robustness of the NFκB response to changes in the feedback parameters (Figure 5). We note that although previous analyses have identified the sensitivity of the NFκB response to many of the same parameters identified here [30,32], none appear to have interpreted the importance of such parameters in the context of feedback control systems. The behavior of the NFκB regulatory network is not unlike that commonly encountered in feedback systems in the engineering world. Consider, for instance, the operation of an amplifier designed to amplify signals in an electronic system. High gain amplifiers with negative feedback amplify signals robustly even when subjected to relatively large changes in feedforward system parameters. But the response is sensitive to feedback parameters, which both permits the system to be finely tuned by selecting proper feedback components, and makes the system vulnerable to failure if the feedback parameters are altered significantly, perhaps as a result of severe damage. This exemplifies the "robust yet fragile" response that is a general characteristic of complex systems with feedback regulation, whether in biology or engineering [46].
In the NFκB signaling network, feedback from IκBαinduced transcription allows the system to respond robustly to stimuli to control gene expression, but at the same time makes the system sensitive to changes in feedback parameters. The highly responsive nature of the system makes it particularly susceptible to network perturbations affecting the feedback molecules IκBα and A20, perhaps as might be seen with severe injury such as stroke. However this feature also provides great opportunities for targeted treatment or intervention to modulate the response. Mathematical modeling and analysis may prove indispensible for future exploration of the NFκB response and drug targeting in microglia, especially when considering crosstalk among multiple pathways that are simultaneously activated by brain injury.
Conclusions
Mathematical modeling has been used extensively in recent years to provide a detailed view into the activation of NFκB, helping to make sense of the multiple layers of feedback and to provide a much deeper understanding of how the system functions as a whole. Here we present the development of a mathematical model that quantitatively describes canonical IKK and NFκB activation in a novel cell type: microglia. The approach we used in model development exploits the multiple feedback structure of the network, and allows the model to be developed in multiple stages by breaking individual feedback loops and developing the modules using the appropriate experimental data. This approach may also prove useful for modeling other biological systems with feedback regulation.
This mathematical model differs significantly from existing NFκB signaling models in two regards. First, it introduces nonlinearities into the activation and inactivation rates for IKK, which are necessary to reproduce the quantitative IKK profile obtained experimentally and correspond with known biological mechanisms. Secondly, the model includes intermediate dynamics in the induced IκBα degradation pathway. We showed these additional dynamics are essential to characterize NFκB signaling observed in microglia in a statistically significant manner and are likely due to reactions involved in the ubiquitination and proteasomal degradation of IκBα, suggesting a more prominent role for this system in modulating the NFκB response.
The mathematical model developed here is the first of its kind for microglia and offers a valuable new tool to study inflammatory signaling in this cell type, permitting rapid numerical simulation and analysis. Our numerical analyses emphasize the highly dynamic nature of regulation of the NFκB network in response to TNFα stimulus, an aspect which has received relatively little attention in prior analyses. While several key parameters play a significant role in modulating the response throughout the entire duration, many others only regulate the response during specific time intervals, such as during the initial activation phase or the oscillatory later phase. The analysis further provides insight into the robustness properties of the system, indicating high sensitivity to feedback parameters, which we note is analogous to the operation of negative feedback systems in engineering.
Methods
Cell culture
BV2 cells, a mouse microglia cell line and kind gift from Dr. K. Andreasson at Stanford University, were cultured in Dulbecco's Modification of Eagle's medium (DMEM, GIBCO, cat# 11995) supplemented with 8% Fetal Bovine Serum (Hyclone, cat#SV30014.03), Penicillin (100 U/ml, GIBCO, cat#15140), and Streptomycin (100 μg/ml, GIBCO, cat#15140). Cells were passaged every four days and were used between passages 1020.
Measurement of activated NFκB p65
BV2 cells were seeded at 4 × 10^{5 }cells per well in six well plates 36 hrs prior to treatment with 10 ng/ml recombinant mouse TNFα (R&D systems, cat# 410MT). Cells were then harvested for protein at the indicated time points with Phosphosafe Extraction buffer (Novagen, cat#71296) supplemented with 0.01 volume Protease Inhibitor cocktail (Sigma, cat# p8340) and 5 mM DTT before use. Protein concentration was measured using the Coomassie Plus assay (Pierce, cat#23236). 25 μg total protein from each sample was transferred to a prechilled Eppendorf tube and brought to 25 μl with complete lysis buffer. These aliquots were stored at 80°C until use for activated NFκB p65 measurement. Active NFκB was measured using the Trans AM NFκB p65 Transcription Factor Assay Kit (Active Motif, cat#40096) according to the manufacturer's instructions. 20 μg total protein was used for each sample. Three cultures were assayed for each group. Standards were prepared from recombinant p65 (Active Motif, cat#31102).
IKK measurements
IKK activity was measured by immunoprecipitation of IKK trimers, followed by a kinase assay/ELISA using a modification of the KLISA IKKβ Inhibitor Screening Kit (Calbiochem, cat# is CBA044). A total of 400 μg protein from each sample was incubated at 4°C 5 hrs with 5 μg goat antiIKKγ antibody M18 (Santa Cruz Biotechnology, Cat# is SC8256) with shaking, followed by overnight incubation with shaking with 50 μl 2 × diluted Protein GSepharose (Sigma, Cat# is P3296) previously washed in complete lysis buffer. Beads were then centrifuged for 5 min at 13,000 rpm 4°C, the postimmunoprecipitation supernatant removed, and beads were washed in the 1 × kinase assay buffer from the KLISA kit. Beads were then incubated with shaking in an incubator for 1 h at 30°C in 75 μl 1 × kinase assay buffer containing 150 ng GSTIκBα and 1 × ATP/MgCl_{2 }mix from the kit. Beads were then centrifuged at 13,000 rpm for 5 min at 4°C, and 60 μl of supernatant was transferred to a well of the glutathione coated 96well plate provided with the KLISA kit. Twofold serial dilutions of the recombinant IKKβ provided with the kit were run as standards according to the kit instructions, but omitting IKK inhibitor. In addition the postimmunoprecipitation supernatant was concentrated 20 × and run to demonstrate that all IKK activity was depleted from the supernatant. In all cases this sample showed no IKK activity. The plate was incubated 30 min at 30°C to allow the GSTIκBα to bind, and subsequent processing was done according to the vendor's instructions. Final concentrations measured were normalized to the total amount of protein used in a given experiment.
Total IκBα measurement
Total IκBα measurements from TNFα treated BV2 cells were performed using the PathScan Total IκBα Sandwich ELISA kit from Cell Signaling (#7360). BV2 cells from passage 1418 were seeded at 4 × 10^{5 }cells/ml on day one and treated with 10 ng/ml TNFα on day three. Cell lysates were prepared and ELISA analysis performed following the manufacturer's instructions. Total protein concentrations were measured using the BCA method; 275 μg total protein was used to measure total IκBα at each time point. The experiments were repeated 3 times.
Analysis of experimental data
Data from each experiment for NFκB and IKK was normalized relative to the maximum mean level of activity during that particular experiment to account for variations in optical absorbance readings between experiments. The normalized data were then averaged to produce the ensemble average data set used for data fitting.
Mathematical modeling and simulation
The model, based on the ordinary differential equation twofeedback model in [14], was developed to incorporate intermediate steps involved in the ubiquitination and proteasomal degradation of IκBα, A20 feedback at multiple points, and nonlinear IKK activation and inactivation rates. The model was integrated numerically using MATLAB 7.7.0 (MathWorks) following the simulation protocol used in [14]. Briefly, the system was initialized with concentrations of total NFκB and IKK, with all other species set to zero (Additional file 1: Table S1). The model was simulated without stimulus for sufficient time to equilibrate the system. Equilibrium concentrations were then used as the initial conditions for simulations with TNFα stimulus present. Active IKK was assumed to be zero during equilibration and to remain constant at a low level of activity at time points beyond 30 min [20,22] for simulations in which the experimental IKK curve was used as input. The IKKa concentration was computed at each time point during simulation using piecewise cubic Hermite interpolation ('pchip') with the interp1 function in Matlab. Similarly, nuclear NFκB was interpolated in an identical procedure from a simulated curve for development of the upstream module. Further details about the mathematical modeling and tables listing all model species, reactions and parameters can be found in Additional file 1 and Additional file 2. The Matlab source code for the ODE model and simulation script are available upon request.
Additional file 2. SBML model of microglial NFκB activation. This .xml file is an SBML translation of the mathematical model originally developed in Matlab (MathWorks). Simulations correspond to the response of microglia cells following treatment with 10 ng/ml TNFα.
Format: XML Size: 30KB Download file
Statistical evaluation of model simulations
The agreement between model simulations and experimental data was assessed using an approach based on Fisher's combined probability test [47], which is justified as follows. Each experimental sample is assumed to be the sum of the population mean and measurement noise. The measurement noise is assumed to be iid Gaussian with zero mean. Each data point itself is the bulk average of a large number of cells (>10^{5}), and so it is assumed that the sample average from this large collection of cells is normally distributed with mean equal to the population average, but that the standard deviation can vary with time. Individual samples are assumed to be independent across experiment replicates and identically distributed with regard to their respective time points. This is justified since all samples are collected from independent cell populations.
Under these assumptions, a twosided one sample ttest can be used to compare the population mean from the model simulations corresponding to a specific set of parameters, θ, to the sample mean from n_{i }experimental samples collected at time t_{i}. The null hypothesis that the two are consistent is rejected at a significance level α if the pvalue corresponding to the ith tstatistic is p_{i}< α.
Fisher's method combines the information from the individual test results to test the shared null hypothesis that all the n_{i }experimental samples come from cell populations whose time evolution of the population average is given by the kinetic model. The test statistic for Fisher's method is computed by combining each independent test as follows:
where log denotes the natural logarithm, and p_{i }are the pvalues from the ttests at n time points. Under the null hypothesis, χ^{2}_{F }follows a chisquare distribution with 2n degrees of freedom. The shared null hypothesis is rejected at a significance level α_{F }confidence if p_{F}< α, thus giving a statistical basis upon which a candidate parameter set can be rejected or retained.
Parameter estimation and sensitivity analysis
Parameter values were estimated by minimizing the a cost function based on the goodness of fit between model and data. Two objective functions were used: one which computed the normalized sum of squares error (SSE),
between the model simulations at parameter set θ, y(t_{i},θ), and observed data points y_{obs}(t_{i}), where i indexes the n time points at which data was collected. A second objective function used the chisquare test statistic computed from Fisher's method (χ^{2}_{F}), an adaptation of the momentmatching algorithm proposed in [23]. The simulated concentrations of NFκB and IKK were normalized to their respective concentrations at 20 min and 5 min to allow direct comparison with experimental data. Optimization was performed using the fmincon constrained minimization algorithm from the Matlab Optimization Toolbox (MathWorks). Lower and upper bounds for the parameter values were taken from the available literature, as specified in Additional file 1.
The normalized first order sensitivity coefficients of the system,
where y_{i }is a system output and θ_{j }is the jth rate parameter, were solved using the CVODES forward sensitivity solver from the SUNDIALS 2.4.0 software suite (Lawrence Berkeley National Labs). Sensitivity scores were also assigned based on the timeaveraged integral of the normalized sensitivity magnitudes,
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
MK, RGG and PWS contributed to the design of experiments, PWS developed the mathematical model and performed the simulations, JFE and XS perfomed the experiments with BV2 cells, PWS, RGG and MK contributed to writing the manuscript. All authors read and approved the final manuscript.
Acknowledgements
We thank Dr. Katrin Andreasson for BV2 cells and Gabriele Lillacci for assistance with the statistical methods and for reviewing the manuscript. This work was supported by NIHR01 GM04983 to RGG and MK and by an NSF IGERT Traineeship DGE0221715 to PWS.
References

O'Neill LA, Kaltschmidt C: NFκB: a crucial transcription factor for glial and neuronal cell function.
Trends Neurosci 1997, 20(6):252258. PubMed Abstract  Publisher Full Text

Schneider A, MartinVillalba A, Weih F, Vogel J, Wirth T, Schwaninger M: NFkB is activated and promotes cell death in focal cerebral ischemia.
Nat Med 1999, 5:554559. PubMed Abstract  Publisher Full Text

John GR, Lee SC, Brosnan CF: Cytokines: powerful regulators of glial cell activation.
Neuroscientist 2003, 9(1):1022. PubMed Abstract  Publisher Full Text

Hoffman A, Levchenko A, Scott ML, Baltimore D: The IκBNFκB Signaling Module: Temporal Control and Selective Gene Activation.
Science 2002, 298:12411245. PubMed Abstract  Publisher Full Text

Nelson D, Ihekwaba A, Elliott M, Johnson J, Gibney C, Foreman B, Nelson G, See V, Horton C, Spiller D, et al.: Oscillations in NFκB Signaling Control the Dynamics of Gene Expression.
Science 2004, 306(5696):704708. PubMed Abstract  Publisher Full Text

Tay S, Hughey JJ, Lee TK, Lipniacki T, Quake SR, Covert MW: Singlecell NFκB dynamics reveal digital activation and analogue information processing.
Nature 2010, 466:267271. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hayden MS, Ghosh S: Shared Principles in NFκB Signaling.
Cell 2008, 132(3):344362. PubMed Abstract  Publisher Full Text

Chen Z, Hagler J, Palombella V, Melandri F, Scherer D, Ballard D, Maniatis T: Signalinduced sitespecific phosphorylation targets I kappa B alpha to the ubiquitinproteasome pathway.
Genes Dev 1995, 9(13):15861597. PubMed Abstract  Publisher Full Text

Suzuki H, Chiba T, Kobayashi M, Takeuchi M, Suzuki T, Ichiyama A, Ikenoue T, Omata M, Furuichi K, Tanaka K: IκBa Ubiquitination Is Catalyzed by an SCFlike Complex Containing Skp1, Cullin1, and Two FBox/WD40Repeat Proteins, βTrCP1 and βTrCP2.
Biochem Biophys Res Commun 1999, 256(1):127132. PubMed Abstract  Publisher Full Text

Yaron A, Hatzubai A, Davis M, Lavon I, Amit S, Manning AM, Andersen JS, Mann M, Mercurio F, BenNeriah Y: Identification of the receptor component of the IκBaubiquitin ligase.
Nature 1998, 396(6711):590594. PubMed Abstract  Publisher Full Text

Alkalay I, Yaron A, Hatzubai A, Orian A, Ciechanover A, BenNeriah Y: Stimulationdependent I kappa B alpha phosphorylation marks the NFkappa B inhibitor for degradation via the ubiquitinproteasome pathway.
Proc Natl Acad Sci USA 1995, 92(23):1059910603. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pahl H: Activators and target genes of Rel/NFkappaB transcription factors.
Oncogene 1999, 18(49):68536866. PubMed Abstract  Publisher Full Text

Wertz IE, O'Rourke KM, Zhou H, Eby M, Aravind L, Seshagiri S, Wu P, Wiesmann C, Baker R, Boone DL, et al.: Deubiquitination and ubiquitin ligase domains of A20 downregulate NFkB signalling.
Nature 2004, 430(7000):694699. PubMed Abstract  Publisher Full Text

Ashall L, Horton CA, Nelson DE, Paszek P, Harper CV, Sillitoe K, Ryan S, Spiller DG, Unitt JF, Broomhead DS, et al.: Pulsatile Stimulation Determines Timing and Specificity of NFκBDependent Transcription.
Science 2009, 324(5924):242246. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Basak S, Kim H, Kearns JD, Tergaonkar V, O'Dea E, Werner SL, Benedict CA, Ware CF, Ghosh G, Verma IM, et al.: A Fourth IκB Protein within the NFκB Signaling Module.
Cell 2007, 128(2):369381. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Covert MW, Leung TH, Gaston JE, Baltimore D: Achieving Stability of LipopolysaccharideInduced NFκB Activation.
Science 2005, 309:18541857. PubMed Abstract  Publisher Full Text

Krishna S, Jensen MH, Sneppen K: Minimal model of spiky oscillations in NFκB signaling.
Proc Natl Acad Sci USA 2006, 103(29):1084010845. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lipniacki T, Paszek P, Brasier AR, Luxon B, Kimmel M: Mathematical model of NFκB regulatory module.
J Theor Biol 2004, 228(2):195215. PubMed Abstract  Publisher Full Text

Cheong R, Hoffmann A, Levchenko A: Understanding NFkB signaling via mathematical modeling.

Cheong R, Bergmann A, Werner SL, Regal J, Hoffmann A, Levchenko A: Transient IκB Kinase Activity Mediates Temporal NFκB Dynamics in Response to a Wide Range of Tumor Necrosis Factorκ Doses.
J Biol Chem 2006, 281(5):29452950. PubMed Abstract  Publisher Full Text

Karin M, BenNeriah Y: Phosphorylation meets ubiquitination: the control of NFκB activity.
Annu Rev Immun 2000, 18(1):621663. PubMed Abstract  Publisher Full Text

Werner SL, Barken D, Hoffmann A: Stimulus Specificity of Gene Expression Programs Determined by Temporal Control of IKK Activity.
Science 2005, 309(5742):18571861. PubMed Abstract  Publisher Full Text

Lillacci G, Khammash M: Parameter estimation and model selection in computational biology.
PLoS Comput Biol 2010, 6(3):e1000696. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Park SG, Lee T, Kang HY, Park K, Cho KH, Jung G: The influence of the signal dynamics of activated form of IKK on NFκB and antiapoptotic gene expressions: A systems biology approach.
FEBS Letters 2006, 580:822830. PubMed Abstract  Publisher Full Text

Werner SL, Kearns JD, Zadorozhnaya V, Lynch C, O'Dea E, Boldin MP, Ma A, Baltimore D, Hoffmann A: Encoding NFκB temporal control in response to TNF: distinct roles for the negative regulators IκBα and A20.
Genes Dev 2008, 22(15):20932101. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lipniacki T, Puszynski K, Paszek P, Brasier AR, Kimmel M: Single TNFα trimers mediating NFκB activation: stochastic robustness of NFκB signaling.
BMC Bioinformatics 2007, 8:376. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Delhase M, Hayakawa M, Chen Y, Karin M: Positive and negative regulation of IκB kinase activity through IKKβ subunit phosphorylation.
Science 1999, 284(5412):309313. PubMed Abstract  Publisher Full Text

Mauro C, Pacifico F, Lavorgna A, Mellone S, Iannetti A, Acquaviva R, Formisano S, Vito P, Leonardi A: ABIN1 binds to NEMO/IKKγ and cooperates with A20 in inhibiting NFκB.
J Biol Chem 2006, 281(27):1848218488. PubMed Abstract  Publisher Full Text

Zhang SQ, Kovalenko A, Cantarella G, Wallach D: Recruitment of the IKK signalosome to the p55 TNF receptor: RIP and A20 bind to NEMO (IKKγ) upon receptor stimulation.
Immunity 2000, 12(3):301311. PubMed Abstract  Publisher Full Text

Ihekwaba A, Broomhead D, Grimley R, Benson N, Kell D: Sensitivity analysis of parameters controlling oscillatory signalling in the NFκB pathway: the roles of IKK and IκBα.
Systems Biol 2004, 1:93103. Publisher Full Text

Ihekwaba AE, Broomhead DS, Grimley R, Benson N, White MR, Kell DB: Synergistic control of oscillations in the NFκB signalling pathway.
IEE Proc Systems Biol 2005, 152(3):153160. Publisher Full Text

Joo J, Plimpton S, Martin S, Swiler L, Faulon JL: Sensitivity Analysis of a Computational Model of the IKKNFkappaBIkappaBalphaA20 Signal Transduction Network.
Ann N Y Acad Sci 2007, 1115:221239. PubMed Abstract  Publisher Full Text

Yue H, Brown M, Knowles J, Wang H, Broomhead DS, Kell DB: Insights into the behaviour of systems biology models from dynamic sensitivity and identifiability analysis: a case study of an NFκB signalling pathway.
Mol Biosyst 2006, 2(12):640649. PubMed Abstract  Publisher Full Text

Lee TK, Denny EM, Sanghvi JC, Gaston JE, Maynard ND, Hughey JJ, Covert MW: A Noisy Paracrine Signal Determines the Cellular NFκB Response to Lipopolysaccharide.
Sci Signal 2009, 2:ra65. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Paszek P, Ryan S, Ashall L, Sillitoe K, Harper CV, Spiller DG, Rand DA, White MRH: Population robustness arising from cellular heterogeneity.
Proc Natl Acad Sci USA 2010, 107(25):1164411649. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Karin M: The beginning of the end: IκB kinase (IKK) and NFκB activation.
J Biol Chem 1999, 274(39):2733927342. PubMed Abstract  Publisher Full Text

Pierce NW, Kleiger G, Shan S, Deshaies RJ: Detection of sequential polyubiquitylation on a millisecond timescale.
Nature 2009, 462(7273):615619. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Nakayama K, Hatakeyama S, Maruyama S, Kikuchi A, Onoé K, Good RA, Nakayama KI: Impaired degradation of inhibitory subunit of NFκB (IκB) and βcatenin as a result of targeted disruption of the βTrCP1 gene.
Proc Natl Acad Sci USA 2003, 100(15):87528757. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Fuchs S, Spiegelman V, Kumar K: The many faces of βTrCP E3 ubiquitin ligases: Reflections in the magic mirror of cancer.
Oncogene 2004, 23(11):20282036. PubMed Abstract  Publisher Full Text

Read MA, Brownell JE, Gladysheva TB, Hottelet M, Parent LA, Coggins MB, Pierce JW, Podust VN, Luo RS, Chau V, alombella VJ: Nedd8 modification of cul1 activates SCF^{βTrCP}dependent ubiquitination of IκBa.
Mol Cell Biol 2000, 20(7):23262333. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kawakami T, Chiba T, Suzuki T, Iwai K, Yamanaka K, Minato N, Suzuki H, Shimbara N, Hidaka Y, Osaka F, et al.: NEDD8 recruits E2ubiquitin to SCF E3 ligase.
EMBO J 2001, 20(15):40034012. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lee EG, Boone DL, Chai S, Libby SL, Chien M, Lodolce JP, Ma A: Failure to Regulate TNFInduced NFκB and Cell Death Responses in A20Deficient Mice.
Science 2000, 289(5488):23502354. PubMed Abstract  Publisher Full Text

SchomerMiller B, Higashimoto T, Lee YK, Zandi E: Regulation of I kappa B Kinase (IKK) Complex by IKK gammadependent Phosphorylation of the Tloop and C Terminus of IKK beta.
J Biol Chem 2006, 281(22):15268. PubMed Abstract  Publisher Full Text

Kearns JD, Basak S, Werner SL, Huang CS, Hoffmann A: IκBε provides negative feedback to control NFκB oscillations, signaling dynamics, and inflammatory gene expression.
J Cell Biol 2006, 173(5):659664. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sillitoe K, Horton C, Spiller D, White M: Singlecell timelapse imaging of the dynamic control of NFκB signalling.
Biochem Soc Trans 2007, 35:263266. PubMed Abstract  Publisher Full Text

Csete ME, Doyle JC: Reverse Engineering of Biological Complexity.
Science 2002, 295(5560):16641669. PubMed Abstract  Publisher Full Text

Fisher RA: Statistical methods for research workers. 11th edition. Edinburgh: Oliver and Boyd; 1950.