Dual Delayed Feedback Provides Sensitivity and Robustness to the NF-κB Signaling Module

Many cellular stress-responsive signaling systems exhibit highly dynamic behavior with oscillatory features mediated by delayed negative feedback loops. What remains unclear is whether oscillatory behavior is the basis for a signaling code based on frequency modulation (FM) or whether the negative feedback control modules have evolved to fulfill other functional requirements. Here, we use experimentally calibrated computational models to interrogate the negative feedback loops that regulate the dynamic activity of the transcription factor NF-B. Linear stability analysis of the model shows that oscillatory frequency is a hard-wired feature of the primary negative feedback loop and not a function of the stimulus, thus arguing against an FM signaling code. Instead, our modeling studies suggest that the two feedback loops may be tuned to provide for rapid activation and inactivation capabilities for transient input signals of a wide range of durations; by minimizing late phase oscillations response durations may be fine-tuned in a graded rather than quantized manner. Further, in the presence of molecular noise the dual delayed negative feedback system minimizes stochastic excursions of the output to produce a robust NF-B response.


Introduction
Many important signal transduction pathways contain a negative feedback motif consisting of an activator that activates its own repressor.Activated repression is capable of generating oscillatory behavior [1] and has been observed to do so in biological systems such as the Hes1 regulatory protein which controls neuronal differentiation [2], the p53-Mdm2 system that mediates the DNA damage response [3], and the NF-kB (Q04207) signaling network that governs the immune response and inflammation [4,5].
The role of activated repression is well understood in the context of transient signaling as functioning to limit the duration of the induced activity.Indeed, misregulation of the negative feedback mechanisms that control NF-kB and p53 has been shown to generate prolonged inflammatory or genotoxic stress responses, respectively, that lead to cell death or chronic disease [6,7].Further, negative feedback can sensitize and speed-up responses to weak or transient input signals [8] when compared to constitutive attenuation mechanisms.
In contrast, the physiological role of oscillatory signaling behavior remains poorly understood.Recent work has shown that, in the calcium stress pathway in yeast, the frequency of nuclear localization of a stress-response transcription factor can be modulated by the magnitude of the extracellular calcium concentration, and this frequency modulation results in a coordinated expression of target genes [9].In the NF-kB and p53 signaling systems, the function of oscillations is still unknown.Oscillations in p53 activity were proposed to represent a counting mechanism that quantizes the response, ensuring a robust but appropriate amount of activity for a specific degree of DNA damage [10].An alternate view was proposed in which oscillations of the p53-controlling ATM kinase activity allow for periodic sampling of the damaged DNA to track its repair and, if necessary, drive further p53 signaling to sustain the repair programs [11].Oscillations in NF-kB activity were proposed to determine which genes would be transcriptionally induced, thereby representing a temporal code that conveys information about the stimulus to gene promoters [5].However, it is not clear whether or not the frequency encodes information in this systems as no differences in NF-kB target gene expression were observed between oscillating and non-oscillating genetic variants [12].
Recent work has demonstrated that oscillations in NF-kB activity can be generated by pulsatile stimulation with TNFa (P06804) [13].However, an analysis of the repeated activation of NF-kB that is driven by an oscillating signal provides little information about the role of oscillations that naturally arise with persistent stimulation.Thus, the role(s) of oscillations in NF-kB activity remains unclear and several questions are still unanswered: Do these oscillations convey information encoded in the frequency to downstream processes?Do they function to generate a periodically recurring phase of sensitivity to stimuli or regulatory crosstalk representing a potential ''counting'' mechanism?Do they ''quantize'' the output signal, thus specifying robust units of activity?Or, are the oscillations caused by persistent signaling simply a non-functional by-product of the requirement for the negative feedback architecture to enable sensitive, fast responses to transient stimuli?
Mathematical models comprised of a small number of equations have led to a greater understanding of biological processes in terms of molecular interactions, diffusion, dose responses, gradient sensing, the role stochasticity in gene expression and in fate decisions [14][15][16][17].Although several models of networks with autoregulation have been developed [18][19][20], most of these networks do not incorporate delays.In signaling, however, such elegant models often do not faithfully reproduce the dynamic behavior of the signaling system because actual biological networks involve many molecular interactions that tend to slow overall signal processing.Larger models comprised of many molecular species and parameters have proven useful in exploring dynamic signaling behavior via computational simulations in conjunction with experimental studies, but they are analytically intractable and therefore do not provide the degree of conceptual insights that small models do.
Here we pursue an alternative approach to modeling NF-kB signaling.We construct a new model that replaces cascading reactions with a single but delayed compound reaction that enables both recapitulation of experimentally observed dynamics and the use of powerful analytical tools.With these tools, we explore the physiological function of the dynamic behavior of NF-kB produced by the activated repression mechanism mediated by its inducible inhibitors, IkBa (Q9Z1E3) and IkBE (O54910).The mathematical analysis results in predictions that are addressed experimentally and thus lead to fundamental insights about the function and origins of this signaling system.

NF-kB model formulation
The basic structure of the NF-kB signaling module is shown in Figure 1 A [4].In resting cells, NF-kB is sequestered in the cytoplasm by IkB proteins.Cellular stimulation leads to activation of the IkB kinase (IKK) which phosphorylates IkB proteins thus targeting them for degradation.Upon degradation of IkB proteins, NF-kB moves into the nucleus and activates hundreds of target genes including the predominant IkB isoform, IkBa.Synthesized IkBa enters the nucleus, binds to NF-kB, and the IkBa-NF-kB complex is exported back to the cytoplasm.Thus, the core feature of the NF-kB signaling module is a negative feedback loop mediated by IkBa.This can be abstracted to a simple motif in which x (NF-kB) activates y (IkBa), y represses x, and repression of x by y is relieved by K (active IKK) (Figure 1B).
Using this motif as a guide, we formulated our model of the IkBa-mediated NF-kB response as a set of 9 reactions and 6 variables (Tables 1, 2).Specifically, the model assumes that the total number of the NF-kB molecules (X ) is conserved, however they can exist either in free/nuclear form (x) or sequestered outside of nucleus within the IkBa-NF-kB complex (½xy).The model contains non-delayed reactions for the binding of free NF-kB to the unbound IkBa promoter (d 0y ) to form the bound IkBa promoter (d 1y ), binding of IkBa protein (y) to free NF-kB to form the IkBa-NF-kB complex, constitutive degradation of IkBa, and

Author Summary
Many signaling events are controlled by negative feedback circuits: as a result they are highly dynamic and in some cases show oscillations The presence of oscillations has led to the hypothesis that signaling pathways convey information about the stimulus via the frequency of oscillations and spikes of activity, analogous to frequency modulated (FM) radio signals.One such signaling protein is NF-kB which controls the inflammatory and immune response to cytokines and pathogens.We show here that the topology of the negative circuit does not allow for frequency modulation by the signaling input.Instead, we show that a second negative feedback circuit may be tuned to dampen the oscillations.In fact, the resulting dual negative feedback motif allows for better tracking of the duration of the incoming signal than the single negative feedback circuit, as well as better buffering of noise present in the incoming signal.Thus we propose that the negative feedback topology has evolved to provide complex dynamics of NF-kB in vertebrate animals and not for the purposes of oscillations.induced degradation of free and bound IkBa proteins by the active IkB kinase IKK (K) producing free NF-kB.In contrast, a compound delayed reaction describes the synthesis of IkBa protein.This reaction involves a time delay t y , which represents the time needed for transcription, translation, nuclear import and export, and protein-protein interactions.
Using experimentally validated assumptions, we reduced the set of mass-action kinetics equations for the 9 reactions to a single delay-differential equation: where Y (y)~y½1zc y X =(1zc y y) is the total IkBa concentration (the sum of free IkBa (y) and IkBa bound to NF-kB), d 0y (y)~Fd y =½F zX =(1zc y y), d 1y (y)~Xd y =½F (1zc y y)zX , are the probabilities for the IkBa promoter to be free or bound to NF-kB, respectively, F~f 1 =f 0 , c y ~ky =(k {y zKr y ), and the subscript t denotes the variable taken at time t{t y (see Methods for details of the derivation).The rates of individual reactions k y ,k {y ,f 1 ,f 0 ,b y ,g y ,r y ,a y are defined in Table 2.
Mirroring the biological system, the non-dimensional timedependent parameter K(t), which characterizes the active IKK concentration, is used as the proxy input signal.The first term in the r.h.s. of Eq. 1 represents constitutive synthesis from the unbound IkBa promoter, the second term represents induced synthesis from the NF-kB-bound IkBa promoter, the third term represents constitutive degradation of IkBa protein, and the fourth term represents IKK-induced degradation of IkBa.Values of K*1 correspond to the rate of IKK-induced degradation of NF-kB-IkBa complex which is of the same magnitude as unbound IkBa.Nuclear NF-kB level x at any time can be determined directly from IkBa levels via x~X =(1zc y y).The time delay t y is incorporated in the synthesis terms: we assume that the rate of production of new proteins at time t depends on the state of the system at time t{t y .Incorporating this time delay allows us to explore the behavior of the negative feedback loop without simulating the full set of reactions associated with it.We obtained values for the time delay and for the other model parameters by calibrating the behavior of the model with experimental results (Table S1).As a starting point, we used parameter values from biochemical measurements [21].However, some modifications were necessary because these values represent the rates of single reaction steps and the model contains compound reactions.
To validate the model, we compared it to experiments.In response to a persistent input signal (starting at time t~0), our simulations of the IkBa-mediated negative feedback system show pronounced oscillations in nuclear NF-kB levels with an oscillation period of about 90 minutes (Figure 1 C).Oscillations with a similar period were observed experimentally when mutant cells containing only the IkBa feedback loop were persistently stimulated with the inflammatory cytokine TNF (Figure 1D).
To address the dynamics of the wild-type NF-kB system that feature both IkBa and IkBE feedback loops, we expanded the model to include an additional 9 reactions and 4 variables involving IkBE (Tables 3, 4).Following the same reduction procedure (see Methods for derivation), we derived a deterministic model consisting of two coupled delay-differential equations for the concentrations of the two IkB isoforms, IkBa (y) and IkBE (z),    S1).
Is the oscillation period a function of the stimulus?
The advantage of our modeling approach is that it allows for analytical studies of the network dynamics.Here, we perform a linear stability analysis of the delay-differential equation (1) to identify the characteristic period and decay rate of NF-kB oscillations produced when input signal is present (Kw0).For sufficiently large K, induced synthesis and degradation are much stronger than basal ones, so the latter can be neglected (a y ~gy ~0).
Expressing y via Y and substituting it into d 0y , d 1y yields a closed equation for Y in the form where Y t ~Y (t{t y ) and the function G( : ) has the form The fixed point Y s (stationary solution) of this equation is given by the algebraic equation The stability of this solution is determined by the eigenvalue of the linearized equation ( The imaginary part of l gives the oscillation frequency v~2pf , and the (negative) real part of l gives the decay rate d of oscillations.Plotting the period (T~1=f ) (Figure 2A) and decay (d) (Figure 2B) of the oscillations as a function of the delay reveals a strong dependence.In contrast, the signaling perturbation K (the active IKK kinase) that acts as the input for the model determines the amplitude of the response but only negligibly affects the period or the oscillation decay (Figure 2B).The mathematical reason for this asymmetry is that the imaginary part of the Lambert function W (x) for negative values of its argument changes very weakly for arguments below {2 (ImW ({2)~1:67:::, ImW ({20)~2:27:::) and asymptotically approaches p for very large negative values of the argument.This is why the period of dampened oscillations (2pt y =Im½W (Bt y e Kryty )) depends strongly on delay t y and only very weakly on K.Meanwhile, the real part of the eigenvalue l (the decay rate) is linearly proportional to K because of the second term in Eq.( 7) and also strongly depends on t y because of the first term.Thus, we find that the period is highly dependent on the delay but is rather insensitive to changes in the input level.This is confirmed by direct simulations of the full nonlinear equation ( 1), where time series of x are plotted for several different values of t y and K (Figure S1).Since variations of stimulus do not lead to significant frequency modulation of NF-kB activity, oscillations of NF-kB are unlikely to encode information about the stimulus.

Damping of oscillations in a dual delayed feedback loop system
The main qualitative difference between the one-loop system considered in the previous section, and the wild-type NF-kB module is the presence of another IkB isoform, IkBE, which also provides negative feedback regulation on NF-kB activity (Figure 3A, B), however with slower kinetics [21].Experimental and computational work has shown that IkBE-mediated feedback can cause damping of IkBa -mediated oscillations [21] and (Figure 3C).More recent computational work has predicted that IkBE-mediated feedback desynchronizes oscillations but does not dampen oscillations in single cells [13].Thus, the mechanism by which IkBE-mediated feedback produces damped oscillations at the population level is not well established.Furthermore, it is unknown whether the damping function of the IkBE-mediated feedback loop has evolved to achieve a specific regulatory function or may simply be a secondary consequence of another function.
We hypothesize that the primary role of the second feedback loop is to mitigate oscillatory behavior produced by the first feedback loop.
To address our hypothesis that IkBE-mediated feedback specifically evolved to dampen IkBa-mediated oscillations, we performed a parameter optimization procedure on the wildtype model (Eqs. 2 and 3) to determine the IkBE synthesis parameters that result in maximum damping.To characterize the degree of damping, we chose the maximum peak-trough difference after 6 hrs as a metric for the persistence of oscillations.According to the definition of this performance metric, ''optimal damping'' occurs when this metric is minimized.In our optimization procedure, we varied two important parameters, the time delay of the second feedback loop t z and the scaling factor m which simultaneously varies the rates of constitutive and induced synthesis of IkBE.Choosing m~0 is equivalent to the complete removal of the IkBE-mediated negative feedback loop while m~1 represents the case in which the inducible synthesis rates for IkBE are the same as for IkBa.The two-dimensional optimization search is shown in a color map (Figure 3D) indicating that the performance metric is minimized at m~0:3,t z ~72.Time course simulations with the optimized parameter set show a high degree of damping (Figure 3G) similar to what is observed experimentally (Figure 3C).To determine whether these optimized parameter values correspond to observations, we measured relevant parameter values experimentally.The synthesis delays for IkBa and IkBE were determined by measuring IkBa and IkBE mRNA levels in a time course of TNF-treated murine embryonic fibroblasts (MEFs) in multiple independent experiments (Figure 3E, S2A,B).The measured delay for IkBa was 25:8+5:4min, and 59:4+12:8min for IkBE, which agrees well with the model prediction for optimal damping.
Since it is difficult to measure the promoter strength experimentally, we employed an implicit way of comparing experiment with the model.To relate the parameter value m to experimental measurements, we set m~0:3 in the model and calculated the ratio of peak values for IkBa and IkBE proteins R m , which we found to be equal 3.9.Then we measured the ratios of basal (unstimulated) to peak protein levels for IkBa and IkBE in experiment via quantitative Western blots of whole cell lysates generated during a TNF time course.These were compared to recombinant protein standards to derive absolute molecule number per cell.Peak IkBa protein levels were measured to be 379,800 molecules per cell, and IkBE 71,300 molecules per cell, with both values being subject to an estimated 25% error (Figure 3F, S2C,D).These protein levels correspond to the experimental peak values ratio R e ~5:3 which is close to the model prediction R m ~3:9.

Duration encoding in a dual delayed negative feedback loop system
We next addressed why the NF-kB signaling module may have evolved to produce oscillatory behavior if the oscillation frequency is not a function of the stimulus and does not constitute a signaling code.We first simulated persistent stimulation of a variant NF-kB system without feedback (we assume that IkBa is constitutively produced, so d 0y ~1, d 1y ~0 in Eq. 1) and found that this system produces long term, nonoscillatory NF-kB activity (Figure 4A Top, blue line).As TNF is secreted in bursts and therefore perceived by surrounding cells as transient or pulse stimulation, we then performed stimulations of pulses 15, 30, and 45 min in duration.In the negative feedbackdeficient NF-kB system, the pulses resulted in transient responses that were attenuated very slowly.Faster attenuation can be achieved by increasing the constitutive synthesis rate, a y .Increasing a y by two orders of magnitude results in pulse NF-kB responses to transient stimuli, but the responsiveness (in amplitude) is much reduced (Figure S4).
We then performed similar simulations in a single negative feedback loop NF-kB system and found that this network topology allows for a rapid shutdown of NF-kB activity for transient inputs (Figure 4A Middle).This suggests that the NF-kB network may have evolved from a pathway without feedback to a pathway with a single negative feedback loop to allow for a more sensitive transient response.Although the negative feedback indeed allows for greater sensitivity, a secondary consequence is that pronounced oscillations arise when the input signal persists for a long time period (Figure 4A Middle, blue line).The addition of a second negative feedback loop with a different time delay can help to dampen these oscillations, while preserving the responsiveness of the signaling module to transient stimuli (Figure 4A Bottom).
By plotting the duration of the response (above a given threshold) we investigated what may be called ''temporal dose response curves'' of the single and dual feedback systems (Figure 4B).The dual feedback system has a response duration close to 60 min for short pulses (v100 min), and a duration proportional to the input duration for longer pulses.The single feedback system has the same behavior as the dual feedback system for short inputs, but for longer inputs the single feedback system produces a quantized response with the same output duration for several different input durations.Our analysis indicates that a dual feedback system is able to produce temporally graded responses, whereas a single feedback system that oscillates does not.Given that the duration of the second phase of the NF-kB response to TNF is a critical determinant of gene expression programs [4], we suggest that the NF-kB system has evolved a dual feedback system that allows for NF-kB activity whose duration is more closely related to the duration of the cytokine stimulus.This fine temporal control, achieved via dual negative feedback, may be critical for complex cytokine-mediated cell-to-cell interactions involved in the adaptive immune response present in vertebrates, but may not be necessary for innate patogen-induced immune responses.We hypothesized that, on an evolutionary timescale, the appearance of dual negative feedback loops that regulate NF-kB activity may coincide with the transition from an innate to an adaptive immune system.To address this hypothesis, we used BLASTP with an E-value cutoff of 1e-25 to search for homologs of the mouse IkBa and IkBE protein sequences in other organisms (see Methods).We found homologs for both IkBa and IkBE, not only in other mammals (such as chimp, dog, platypus), but also in other vertebrate classes including fish, amphibians, and birds (Figure 5).Thus, dual negative feedback regulation of NF-kB activity appears to be present in all organisms with adaptive immunity.In contrast, we did not find any invertebrate organisms with homologs for both IkBa and IkBE (Figure 5).Therefore, invertebrates, which lack adaptive immunity, also appear to lack the potential for dual negative feedback regulation of NF-kB mediated by IkBa and IkBE suggesting that the temporal control achieved with this regulatory architecture is not necessary for innate immune responses.

Robustness to fluctuations in a dual delayed negative feedback loop system
Thus far, we have examined the response of the network to transient stimulation in the absence of fluctuations.However, it is well known that noise in gene expression can cause significant variability in cellular responses [18,[22][23][24][25][26].Sometimes this variability can be beneficial [27], but in most cases, noise has a detrimental effect on the robustness of cellular functions.Mechanisms have presumably evolved to mitigate the unwanted effects of noise, especially in signaling pathways.In this section we examine the variability in the response of the NF-kB module that arises due to intrinsic and extrinsic noise, and we demonstrate that the dual-feedback loop architecture allows for a more robust response than the single feedback loop system.Further, we investigate how the relative contribution of intrinsic and extrinsic fluctuations depends on the size of the system.
The concentration of signaling molecules such as NF-kB can vary significantly between cells [28].This variability in protein levels represents a source of extrinsic noise.We examined the variability in the response of the network to fluctuations in the total level of NF-kB and fluctuations in the IKK input level by simulating the network behavior with total NF-kB levels and active IKK levels distributed within a certain rage around their nominal values.The coefficient of variation (CV) in peak nuclear NF-kB levels and the CV in late-phase nuclear NF-kB levels is defined as CV ~2(x max {x min )=(x max zx min ) where x max (x min ) are the maximum (minimum) values of NF-kB at the peak or during the late phase.NF-kB late-phase response is defined as the nuclear NF-kB level following the trough after the first peak response.In Text S1 we compare the extrinsic CV in the peak and the late phase for various values of IKK and NF-kB (see Figure S5).
Intrinsic noise arises from the stochastic nature of biochemical processes such as transcription and translation [24].To examine the response of the NF-kB signaling module in the presence of intrinsic genetic noise, we used the Gillespie algorithm [29] modified according to [30] to perform stochastic simulations of both regular and delayed biochemical reactions included in our delayed feedback model.These latter reactions are initiated at times dictated by their respective rates, but the numbers of molecules are only updated after the time delay since the reaction initiation.
We ran stochastic simulations of both a single and dual feedback system and estimated the ensemble average SX T of the number of NF-kB molecules X and the magnitude of fluctuations as characterized by the standard deviation DX ~SX {SX TT 2 Â Ã 1=2 and the coefficient of variation CV ~DX =SX T. To determine how the variability in the response varies with the magnitude of the input and the size of the system, we determined the CV in peak nuclear NF-kB levels and the CV in late-phase nuclear NF-kB levels for several values of IKK (Figure 6A,C) and for systems with up to 100,000 NF-kB molecules (Figure 6B,D).In Figure 6, we also plot CV values for extrinsic variations (+20%) in total NF-kB at several values of IKK (Figure 6A,C) and CV values for extrinsic variations in IKK (+20%) for several different system sizes (Figure 6B,D).We find that, even with this relatively low level (+20%) of extrinsic variability in IKK and NF-kB protein levels [28], variability in the response of the network is dominated by extrinsic noise for large systems (w10,000 NF-kB molecules).
The CV in late-phase nuclear NF-kB levels is similar for extrinsic and intrinsic noise when the size of the system is reduced to 1000 NF-kB molecules.Next, we investigated the behavior of the NF-kB signaling module in this regime where intrinsic noise levels become significant by analyzing stochastic simulations produced with a system with total NF-kB levels set to 1000 molecules.We ran stochastic simulations of all three systems studied deterministically above: no-feedback, single negative feedback, and dual negative feedback (Figure 7).Note that ensemble-averaged time series agree with the deterministic simulations very well (Figure S6).In the case of no feedback (Figure 7A) there is a strong robust response to the incoming persistent signal as characterized by the low values of the coefficient of variation.However, as we have seen above in Figure 4A, the major flaw of this system is its slow response to the pulse-like signals.Next, we simulated the 9 biochemical reactions included in the IkBa-mediated single negative feedback loop (Figure 7B).In single runs the first peak in nuclear NF-kB levels appears to be very robust, as illustrated by Figure 7B Top.The CV is lowest (v0:2) during the first peak in nuclear NF-kB indicating that this portion of the response is very robust.Subsequent peaks in this undamped system lead to higher CV (w0:5) in the later portion of the response.
Next, we performed stochastic simulations of the 18 biochemical reactions included in the dual delayed feedback model (with both IkBa-and IkBE-mediated feedback) (Figure 7C).In the dual feedback model, as in the single IkBa-mediated feedback model, there is a very robust first peak.However, unlike the single IkBamediated feedback model, in the dual feedback system the noise levels remain at a low level (v0:5) following the first peak in nuclear NF-kB (Figure 7C Bottom).Thus, the dual feedback architecture allows for lower noise levels also in the later portion of the response.
What is the underlying reason for the robustness of the initial response from this circuit?The main source of intrinsic noise lies in the transcription and translation of IkB isoforms, since they are transcribed from single genes.In contrast, fluctuations in protein degradation and transport processes are relatively small, because the copy numbers of the corresponding molecules are large.In the NF-kB network, the peak in nuclear NF-kB levels that occurs following stimulation is produced via the degradation of IkB proteins that bind and sequester NF-kB in the cytoplasm.Thus, we argue that robustness of the initial response of the NF-kB circuit is explained by the fact that it uses the sequestering mechanism and does not rely on the protein production.
To test this hypothesis, we simulated the behavior of an alternative network that relies on transcription of auto-repressor, rather than the degradation of inhibitor proteins, for signaling (Figure 7D).This system can be modeled with two variables: x, the number of repressor molecules, and d, the binary state of the promoter (d~0 corresponds to the unbound promoter and d~1 corresponds to bound promoter), and with four reactions (binding and unbinding of the repressor to the .Schematic of a phylogenetic tree showing organisms in which IkB homologs were found using BLASTP.Organisms with homologs for both IkBa and IkBE are in blue shaded region and organisms with a single homolog are in red shaded region.The branches in the schematic phylogenetic tree are not drawn to scale.(For simplicity, not all organisms with single or dual homologs are shown here.A complete list is provided in Table S6).doi:10.1371/journal.pcbi.1003112.g005promoter, degradation of the repressor, and delayed synthesis of the repressor with rate K(t)(1{d) where K(t) is the external signal (Tables S2, S3).The input signal activates the production of the auto-repressor which after a certain time delay binds to the promoter and terminates further synthesis.Deterministically, this circuit also provides a desired response to a persistent stimulation with a large well-defined first peak.However, stochastic simulations reveal significant differences in the noise performance of this design as compared with the NF-kB circuit (note that the agreement between deterministic and stochastic simulations is less accurate in this case because of the strong promoter fluctuations (Figure S6D).Activation of the autorepressor network is much less robust than the activation of the NF-kB network (cf. Figure 7D and Figures 7B,C).In fact, in the auto-repressor network, the coefficient of variation is highest (w1) during the initial peak (Figure 7D Bottom).These results confirm our conjecture that the sequestering mechanism incorporated in the design of the NF-kB network gives rise to a much more robust activation of NF-kB than alternative networks that rely on transcription for activation and signaling.This finding is in accord with recent work [31] where the sequestering of Cdc20 protein was also implicated in the noise resistance of the spindle assembly checkpoint.
As we mentioned previously, recent computational work has suggested that persistent oscillations are present in wild-type cells with both IkBa-and IkBE-mediated feedback but stochastic variability leads to desynchronization among individual cells and therefore produces damped oscillations at the population level [13,32].Our computational results demonstrate that, although stochastic oscillations are still present in individual cells with both IkBa-and IkBE-mediated feedback (Figure 7C), the oscillatory propensity can be greatly reduced by the second feedback loop in the wild-type NF-kB signaling module.Further, stochastic simulations of the dual-feedback network reveal highly synchronized damped oscillations (Figure S7C) with cellular variations due to intrinsic noise becoming significant only when the system size is drastically reduced (Figure 7C).
To show that our results are not limited to the conceptual NF-kB model introduced above, we simulated the more detailed stochastic NF-kB model formulated in [32], which explicitly incorporates IKKK/IKK signaling cascade and NF-kB shuttling between the nucleus and the cytoplasm (see Methods and Figure S8A).One of the key assumptions made in the model [32] is that the strong stochasticity of the NF-kB response is caused by the slow and stochastic binding/dissociation of NF-kB to the corresponding promoters of IkBa, IkBE, and A20 target genes.The slow rates chosen by the authors for these reactions lead to the high variability of oscillatory dynamics among cells (Figure S8B).However, there is experimental evidence that the binding time of NF-kB may be significantly shorter, at least in certain types of cells.According to Fluorescence Recovery After Photobleaching (FRAP) measurements in HeLa cells [33], the typical time scale of NF-kB binding to the target promoters is on the order of a second rather than minutes, suggesting more rapid equilibration between the NF-kB-bound promoters and the pool of unbound nuclear NF-kB molecules.We found that increasing the binding and dissociation rates by 10 2 :::10 3 times profoundly changes the dynamics of the signaling system.NF-kB trajectories become more regular, suggesting that the behavior of individual cells translates more directly into the behavior of the population (Figure S8C).After adjusting the binding/dissociation rates along with a few other parameters (Table S5), the updated model recapitulated the population response to chronic TNFa stimulation under various genetic conditions (WT, IkBe {={ , and A20 {={ ) (Figure S9) in agreement with earlier experimental results [4,21,34].
To quantify the magnitude of the late oscillatory NF-kB response to a chronic TNFa stimulation, we chose as a metric the average maximum peak-trough difference 5 hrs after initial stimulation.This quantity can be computed in two different ways.The mean single-cell variability can be characterized by the magnitude M s found by computing the maximum peaktrough differences for individual trajectories, and then averaging them over all trajectories: The population-level variability can be characterized by the magnitude M p which is found by first computing an average trajectory and then computing its maximum peak-trough difference: If the stochasticity is small, these two measures are similar, however for strong stochasticity they may differ significantly.
Using these metrics, we first confirmed that for the parameter values adopted by [32], the model shows significant single-cell oscillations both in the IkBE {={ and in the WT, independently of the time delay in the IkBE loop (M s , Figure 8A), but the population-averaged response shows significant oscillation dampening for the time delay around 45 min (M p , Figure 8C).However, for our re-parameterized model with fast binding/dissociation, the stochasticity of individual trajectories is small, and both metrics show similar trend: the amplitude of oscillations in the WT is strongly suppressed at the optimal time delay of 45 min both for the population average (Figure 8B) and the individual cells (Figure 8D), which falls within the margin of error of our experimental results (Figure 3D).

Discussion
In this work we have developed a minimal model of the NF-kB signaling pathway that uses a small number of reactions (some of them compound) thus making it amenable to mathematical analysis.Previously, another simplified model of NF-kB signaling was developed in which a massive overshoot in IkBa resulted in an effective slowing of signaling dynamics [35], and produced spiky oscillations that are not seen in physiological conditions.Our model, which utilizes an explicit time delay, recapitulates experimentally observed signaling behavior.It demonstrates that models with explicit time delays can be useful for investigating the mechanistic basis of the dynamic behavior of signaling pathways.
Using this model, we explored the potential role of NF-kB oscillations which are observed in a variant of the NF-kB signaling module with the secondary negative feedback loop involving IkBE, disabled.In particular, we addressed the question of whether the frequency of these oscillations contains information, as in neurons which sometimes encode information in the frequency of action potentials [36] and in the activation of the transcription factor NF-AT which is responsive to the number of Ca 2z pulses [37].By analyzing the oscillatory response of a system regulated solely by the IkBamediated negative feedback loop, we found that both the frequency and the decay rate of the oscillations produced by this system are highly dependent on the internal parameters of the circuit, but are not sensitive to changes in the input signal levels.This result suggests that the oscillatory frequency does not encode information about the stimulus.Hence, stimulusspecific gene expression is unlikely determined by stimulusspecific frequencies of NF-kB oscillations.If there is a temporal code for stimulus-specific gene expression it is unlikely to involve frequency modulation, but may involve amplitude modulation over time.
When a second feedback regulator, IkBE, is added to the model, the oscillations caused by a persistent stimulation are significantly dampened, in agreement with our earlier findings [21].By performing an optimization procedure, we determined that the specific experimentally observed parameter values for the synthesis delay and peak protein abundance of both IkB isoforms correspond to maximal efficiency of damping.These findings suggest that the second feedback (IkBE) has evolved to produce damping of the oscillatory Figure 7. Stochastic model simulation results for various network architectures (with 1000 total NF-kB molecules).The architectures analyzed are the NF-kB network with no feedback loops (A), only IkBa-mediated negative feedback (B), the NF-kB network with both IkBa-and IkBEmediated negative feedback (C), and an alternative auto-repressive network (D).The top panel in each group shows four typical runs of stochastic simulations for each network, the middle panel shows the mean and standard deviation for 200 runs of each network, and the bottom panel shows the corresponding coefficient of variation.The input signal, K(t), is switched from K~0 to K~K max at t~20 hrs.In A-C, the magnitude of external signal K max ~2, in D, K max ~50.doi:10.1371/journal.pcbi.1003112.g007behavior of the first feedback (IkBa).Furthermore, we demonstrated that this finding is not limited to our simple model, but can be expanded to more complex models.For example, in a recent model by [32] with fast binding/ unbinding rates of NF-kB the secondary IkBE feedback leads to a reduction in NF-kB oscillations in individual cells.However, cell-cell variability and extrinsic noise can further reduce NF-kB oscillations on a population level.
From the evolutionary perspective, we have a peculiar situation in which a signaling module apparently first developed a negative feedback loop that made it prone to oscillations, and then added a secondary loop which mitigated these oscillations.This brings the question, if oscillatory responses are not beneficial to the cell, why has the primary negative feedback appeared in the system in the first place?By comparing transient response of several variants of signaling modules (0-, 1-and 2-feedback loop designs) in the presence of stochastic fluctuations we showed that the primary negative feedback loop involving the release of sequestered NF-kB proteins created a strong, rapid, and robust response to short pulses of active IKK signal.However, for longer signals a single-feedback-loop system exhibits a suboptimal ''temporal dose response behavior'' that leads to a quantized response to signals of different durations.In contrast, the dual feedback network generates response durations that are proportional to the stimulus input durations.Fine-tuning of the response duration may be reflective of a signaling code in which duration of NF-kB activity may be a key determinant of stimulus-specific gene expression program.
Cytokines such as TNFa facilitate adaptive responses at the effector stages [38].The evolution of cytokines is associated with the evolution of an adaptive immune system to allow for coordination of various cell types [39].Unlike pathogen exposure, cytokines are produced during varying amounts of time thereby generating time-varying signals.Our analysis showed that the dual negative feedback module is more capable at distinguishing differences in the duration of incoming signals.This function is important for the transduction of cytokine signals, but not pathogen signals.Our BLASTP analysis indeed demonstrates that the evolution of the dual negative feedback system may correlate with the evolution of adaptive immunity.

Derivation of the deterministic model
Using mass action kinetics, the full set of reactions for the dual feedback loop NF-kB system (Tables 2, 4) can be expressed by the following ODEs: The total number of kB binding sites on each promoter is conserved: We assume that the total amount of NF-kB in the cell X is conserved Since the number of binding sites available for NF-kB protein is small, we can neglect the amount of NF-kB bound to the IkBa and IkBE promoters, so Solving Eq. 22 for x yields: DNA binding reactions are usually fast, so we can assume that they are at quasi-equilibrium at all times, Using Eqs.19 and 20, substituting into Eqs.24 and 25, and solving for d 0y , d 1y , d 0z , d 1z yields: where F ~f1 =f 0 .We also assume quasi-equilibrium for IkB NF-kB binding reactions, Y and Z can in turn be expressed via y and z by: where c y ~ky =(k {y zKr y ) and c z ~kz =(k {z zKr z ).Equations 34-35 combined with definitions Eqs.26-29, 33, 36, and 37 represent a closed system of two delay-differential equations 2, 3 for the dual-feedback NF-kB module.Setting Z~z~0 in these equations leaves us with a single delay-differential equation for the single feedback loop system Eq. 1.

Details of the linear stability analysis
The fixed point Y s of Eq. ( 4) is given by the algebraic equation (6).Unfortunately, Eq. ( 6) does not permit finding Y s in explicit form.However, this calculation can be significantly simplified if the total number of NF-kB proteins is large, so c y X &1, then y can be neglected as compared with total Y .Then x~X {Y , and d 1y ~dy (X {Y )=(F zX {Y ), and expression (5) for G(Y ) simplifies: Now the stationary level of Y can be obtained explicitly The stability of this stationary solution is determined by the linearized equation ( 4) for a small perturbation j near Y s , _ j j~Bj t {Kr y j ð40Þ where j~Y {Y s , subscript t again indicates the delayed value of j taken at time t{t, and B~dG(Y s )=dY .Using formula (38) we obtain where Y s is given by Eq. (39).The eigenvalue l of the linearized equation ( 40) is found by substituting j~j 0 exp(lt), yielding the transcendental equation whose solution is given by Eq. ( 7).

Stochastic model formulation
For the analysis of a full NF-kB system, we adopted the basic structure of the NF-kB model formulated in [32] which in turn was based on the population-level model first proposed in [4].The structure of the model is shown in Figure S8 A. In resting cells, NF-kB is sequestered in the cytoplasm by IkB proteins.In response to TNFa stimulation, IKKK protein becomes active, and activates IKK kinase.IKK phosphorylates IkB proteins targeting them for degradation.Upon degradation of IkB proteins, NF-kB moves into the nucleus and activates hundreds of target genes.In the model, we focus on the dynamics of three genes associated with the negative feedback of the system.Following NF-kB activation, synthesized A20 proteins attenuate TNFa signal by repressing IKKK and IKK transitions into their active states.NF-kB also binds IkBa and IkBE protein promoters, which following translation in the cytoplasm, translocate back into the nucleus and bind free NF-kB sequestering it out of the nucleus.In addition, IkB proteins are directly responsible for NF-kB dissociation from the DNA.
The biological processes in the model were interpreted through stochastic and deterministic representations similar to [32].Nuclear transport, complex formation, synthesis, transcription, and translation were described through a set of ordinary differential equations (Text S1).Regulation of gene activity through NF-kB binding and dissociation from DNA was modeled using stochastic representation.The time-evolution of the system was accomplished through a hybrid simulation algorithm that uses Gillespie algorithm [29] to evaluate the state of stochastic processes and an ODE solver to compute the state of deterministic processes.

Details of the BLASTP search for IkBa and IkBE homologs
We performed two BLASTP searches (using default parameters) to search for IkBa and IkBE homologs.The mouse IkBa protein sequence (gi28386026) was used as the query for the first search.The mouse IkBE protein sequence (gi2739158) was used as the query for the second search.We used an E-value of 1e-25 as a cutoff for both searches.Homologs for IkBa were found in the organisms listed in Table S6, and homologs for IkBE were found in the organisms listed in Table S7.
Note that we selected only unique homologs for both IkBa and IkBE in all vertebrates.We did not find unique IkBE homologue for several vertebrates.We expect that this is due to the fact that complete genomes are not currently available for these organisms.Table S8 lists the genome status (as of 6/1/11) of all organisms for which IkBa or IkBE homologs were found (http://www.ncbi.nlm.nih.gov/genomes/leuks.cgi).

Cell culture experiments
Immortalized murine embryonic fibroblasts [4] were chronically stimulated with 10 ng/mL TNF (Roche) and IkBa and IkBE mRNA and protein levels were monitored by RNase Protection Assay (RPA) and Western Blot, respectively, as previously described [21].RPA results for each time course were quantitated using ImageQuant software (GE Healthcare) and used to determine the time of half-maximal inducibility between basal and peak mRNA levels for IkBa and IkBE (Figure S2 A,B).Western Blot results were also quantitated with ImageQuant software and used to determine the time point of peak expression.The basal abundances of IkBa and IkBE protein were determined via comparison to a standard curve of recombinant IkB protein (R Tsu, JD Kearns, C Lynch, D Vu, K Ngo, S Basak, G Ghosh, A Hoffmann in preparation).The peak abundances of IkBa and IkBE were determined via multiplication of the basal value by the fold inducibility at the peak time point (Figure S2 C,D).Experimental levels of nuclear NF-kB in cells with only the IkBa-mediated negative feedback loop intact and in wild-type cells containing both IkBa-and IkBE-mediated negative feedback were determined by EMSAs in [4].

Figure 1 .
Figure 1.Oscillatory behavior from a system with a single negative feedback loop.(A) Diagram of the IkBa-NF-kB signaling module.(B) Diagram of a system with a single delayed negative feedback loop.(C) Nuclear NF-kB levels (x) in response to persistent stimulation as a function of time produced using our delayed feedback model.(D) Experimental levels of nuclear NF-kB (determined by EMSAs) in cells with only the IkBa-mediated negative feedback loop intact (data from [4]).doi:10.1371/journal.pcbi.1003112.g001

Figure 3 .
Figure 3. Damped oscillations with a dual negative feedback system.(A) Diagram of the dual feedback system with both IkBa-and IkBEmediated negative feedback loops.(B) Diagram of a system with dual delayed negative feedback loops.(C) Experimental levels of nuclear NF-kB (determined by EMSAs) in wild-type cells containing both IkBa-and IkBE-mediated negative feedback (data from [17]).(D) Optimization of the parameters of the second feedback loop m and t z towards maximizing the oscillations damping.The optimization method minimizes peak-minustrough differences six hours after the onset of stimulation, the global minimum occurs at t z ~72 min, m~:3.The black dot indicates the experimentally measured parameter values (t z ~59 min, m~:2).Note that m was not measured directly.The value of m corresponding to the experimentally measured value of R e was determined with the model (Figure S3).(E) Experimental measurements of IkBa and IkBE synthesis delays.(F) Experimental values for peak IkBa and IkBE protein levels.(G) Simulated time course of nuclear NF-kB levels (x) for the single feedback system and for the optimized dual feedback system in response to persistent stimulation with K~2.doi:10.1371/journal.pcbi.1003112.g003

Figure 4 .
Figure 4. Response of the NF-kB signaling module to transient inputs with magnitude K~2.(A) Time series of x for a system with all feedback removed (top), a system with IkBa-mediated negative feedback (middle), and a system with both IkBa-and IkBE-mediated negative feedback (bottom) in response to 15 min (red), 30 min (orange), 45 min (green), and persistent (blue) stimulation.(B) The response duration as a function of the stimulus duration for the single feedback and dual feedback systems.The response duration is the amount of time x exceeds a threshold level of 50 (as indicated by the dashed black lines in the graphs shown in (A).doi:10.1371/journal.pcbi.1003112.g004

Figure 5
Figure 5. Schematic of a phylogenetic tree showing organisms in which IkB homologs were found using BLASTP.Organisms with homologs for both IkBa and IkBE are in blue shaded region and organisms with a single homolog are in red shaded region.The branches in the schematic phylogenetic tree are not drawn to scale.(For simplicity, not all organisms with single or dual homologs are shown here.A complete list is provided in TableS6).doi:10.1371/journal.pcbi.1003112.g005

Figure 6 .
Figure 6.The coefficient of variation (CV) in nuclear NF-kB levels due to extrinsic and intrinsic fluctuations.The CV was calculated for peak (A,B) and late-phase (C,D) nuclear NF-kB levels for both single and dual feedback systems.The CV due to intrinsic fluctuations was determined from at least 50 runs of the stochastic simulations at each value of IKK (A,C) and total NF-kB (B,D).The CV due to extrinsic fluctuations in total NF-kB and IKK levels was determined by varying the total NF-kB level by +20% for each value of IKK (A,C) and by varying IKK by +20% for value of total NF-kB (B,D).doi:10.1371/journal.pcbi.1003112.g006

Figure 8 .
Figure 8.Effect of delay time on damped NF-kB oscillations in the detailed NF-kB model.Magnitudes of single-cell oscillations 5 hours after stimulation (A,B) and of population averaged oscillations (C,D) are shown for the IkBE {={ knockout (red bar) and in the WT for different time delays in the IkBE feedback loop (blue bars) with original (A,C) and adjusted (B,D) parameter values.The optimal time delay of 45 min is shown by the green bar.Each bar represents the average of nuclear NF-kB variation for 500 single cell trajectories.Error bars represent + one standard deviation.doi:10.1371/journal.pcbi.1003112.g008

Figure
Figure S1 Oscillations produced by the IkBa-mediated negative feedback system.(A) K~2 for t y ~20 min, t y ~30 min, and t y ~40 min and with (B) t y ~25 min for K~:5, K~1, and K~2.(EPS) Figure S2 Representative experimental data for IkBa and IkBE synthesis delays and feedback strengths.(A) mRNA synthesis for IkBa and IkBE were measured by RNase Protection

Table 1 .
~ay d 0y (y t ,z t )zb y d 1y (y t ,z t ){g y y{Kr y Y ð2Þ Single feedback model variables.

Table 2 .
Single feedback model reactions.
doi:10.1371/journal.pcbi.1003112.t004Figure 2. Period and decay rate of oscillations produced by the IkBa-mediated negative feedback system.(A) The oscillation period T as a function of K with t y ~25 min (green line) and as a function of t y with K~2 (red dashed line).(B) The oscillation decay rate d as a function of K with t y ~25 min (green line) and as a function of t y with K~2 (red dashed line).doi:10.1371/journal.pcbi.1003112.g002 .1003112.g008_ x x~{f 0 d 0y xzf 1 d 1y {k y yxzk {y ½xy zKr y ½xy{f 0 d 0z xzf 1 d 1z { k z zxzk {z ½xzzKr z ½xz 0z (t{t z )zb z d 1z (t{t z ){k z zxzk {z ½xz{g z z{Kr z z ð16Þ y yx~k {y ½xyzKr y ½xy ð 30Þ k z zx~k {z ½xzzKr z ½xz ð 31Þ Substituting ½xy and ½xz from Eqs. 30 and 31 into Eq.23 yields: x~X {½k y yx=(k {y zKr y ){½k z zx=(k {z zKr z ) ð32Þ Now we can solve Eq. 32 for x x~X 1zk y y=(k {y zKr y )zk z z=(k {z zKr z ) ð33Þ and substitute it in Eqs. 12 and 16.These equations contain both fast and slow terms.However, it is easy to see that rate equations for variables Y ~yz½xy and Z~zz½xz contain only slow terms: _ Y Y ~ay d 0y (t{t y )zb y d 1y (t{t y ){g y y{Kr y Y Text S1 Additional model details.Extrinsic noise in dual negative feedback loop system and details of the full stochastic model. (PDF)