Oscillation dynamics underlie functional switching of NF-κB for B-cell activation

Transcription factor nuclear factor kappa B (NF-κB) shows cooperative switch-like activation followed by prolonged oscillatory nuclear translocation in response to extracellular stimuli. These dynamics are important for activation of the NF-κB transcriptional machinery, however, NF-κB activity regulated by coordinated actions of these dynamics has not been elucidated at the system level. Using a variety of B cells with artificially rewired NF-κB signaling networks, we show that oscillations and switch-like activation of NF-κB can be dissected and that, under some conditions, these two behaviors are separated upon antigen receptor activation. Comprehensive quantitative experiments and mathematical analysis showed that the functional role of switch activation in the NF-κB system is to overcome transient IKK (IκB kinase) activity to amplify nuclear translocation of NF-κB, thereby inducing the prolonged NF-κB oscillatory behavior necessary for target gene expression and B-cell activation.


INTRODUCTION
Transcription factor nuclear factor kappa B (NF-κB) is activated in response to a variety of extracellular stimuli and regulates transcription of multiple genes involved in cell fate decisions. 1 Given its central role in many cellular processes, the regulation of this transcription factor is critical in several human diseases. [2][3][4][5] In the canonical signaling pathway in B lymphocytes, NF-κB is activated by B-cell receptor (BCR) signaling via a complex chain of events involving protein kinase C β (PKCβ), CARD containing MAGUK protein1 (CARMA1, also known as CARD11), transforming growth factor (TGF) β-activated kinase 1 (TAK1, also known as MAP3K7) and IκB kinase β (IKKβ). In unstimulated cells, NF-κB is retained in the cytosol in an inactive state due to its association with inhibitors of NF-κB (IκBs). After stimulation, phosphorylationinduced degradation of the IκBs by IKKβ liberates NF-κB for nuclear translocation, where it regulates transcription of target genes ( Figure 1a, Supplementary Section 1). NF-κB activity is further regulated by multiple feedback loops in the network, resulting in complex dynamical activity profiles. [6][7][8][9][10][11][12][13] NF-κB activity exhibits two characteristic behaviors, oscillations and switch-like activation, regulated by negative and positive feedback, respectively. Prolonged oscillatory nuclear translocation of NF-κB, whose dynamics are important for induction of gene expression, has been documented in a variety of cells stimulated with ligands such as tumor necrosis factor (TNF). 8,11,14,15 NF-κB also shows cooperative switch-like activation dynamics in response to an increased dosage of TNF ligand at the single cell level 11 and in bulk populations of BCR-activated B cells. 13 It is known that interlinked positive feedback loops associated with negative feedback loops are the basis of sustained periodic oscillations in the cell cycle and in circadian rhythms. [16][17][18][19] However, unlike these examples, positive feedback loops identified in the NF-κB signaling network do not induce those types of NF-κB oscillations, but instead induce damped oscillations. These observations suggest that the positive feedback loops embedded in the NF-κB system have different biological functions from other well-studied biological oscillators. In particular, it is unclear whether switch-like activation and oscillations are separate properties emerging from different regulatory mechanisms in the network, or represent two sides of the same coin.
In this report, we expand our previous mathematical model 13 to include the BCR signaling pathway and transcriptional regulation of NF-κB to analyze the dynamic behaviors of NF-κB in B cells stimulated with extracellular surrogate antigen. We show using the model and quantitative experiments followed by validation analysis that NF-κB activity in B cells exhibits oscillations after switch-like activation. Using a variety of B cells with artificially modified NF-κB signaling networks, oscillations and switch-like activation of NF-κB could be dissected and, in some conditions, these two behaviors were separated. We show that switch-like activation caused by positive feedback loops is essential to compensate for the transient nature of upstream IKK activity and to induce amplified accumulation of nuclear NF-κB, which is in turn necessary for inducing subsequent oscillations for target gene expression.

RESULTS
A comprehensive mathematical model for the NF-κB network To investigate the roles and regulatory mechanisms of feedback loops for NF-κB activation dynamics in BCR signaling, we first constructed a comprehensive mathematical model by integrating earlier models for BCR signaling 13 and transcriptional regulation of NF-κB 9 (Figure 1a). Our experimental results (Figure 1b) revealed an oscillatory response reminiscent of that seen with TNF, 6 an observation suggesting that the NF-κB-IκB system may be regulated similarly in both BCR and TNF signaling networks. We then extended the models, incorporating details of the formation of the CARMA1, B cell chronic lymphocytic leukemia 10 (BCL10), and/or mucosa-associated lymphoid tissue (MALT) lymphoma translocation gene 1 (MALT1) (CBM signalosome) complex, which has critical roles in regulating BCR signaling. Adaptor proteins BCL10, MALT1 and CARMA1, which were implicitly considered but not included in our previous model, 13 were newly added in the current model to examine and validate the effects of those molecules on NF-κB dynamics. The current model includes multiple feedback loops. Two positive feedback loops are   operative in a phosphorylation cascade for TAK1 and IKKβ; one is the enhanced TAK1 activation mediated by IKKβ-dependent CARMA1 phosphorylation at serine 578 (S578), 20 and the other is introduced by trans auto-phosphorylation of IKKβ. 21 The former positive feedback loop reactivates and sustains TAK1 activity upon second phosphorylation of CARMA1 at S578 after the first PKCβ-dependent CARMA1 S668 phosphorylation. 13 These positive feedback regulations contribute to the switch-like activation of IKKβ and NF-κB. In addition, there are three transcriptionally inducible negative feedback loops located in the downstream region of the network. Signal-activated NF-κB translocates into the nucleus and initiates transcription of a variety of genes, including those encoding IκBα and IκBε, which force NF-κB into an inactive state in the cytosol. TNF alpha-induced protein 3 (A20, also known as TNFAIP3) is also induced upon NF-κB activation and negatively regulates IKKβ activity. 22,23 Those negative feedback loops are thought to be important for induction of the oscillatory response of NF-κB. Model parameters were optimized to recapitulate experimental observations obtained using chicken DT40 B cells stimulated with anti-IgM (mAb M4; details of the model are described in Supplementary Figures 1-4 Figure 1c). The Hill coefficient, which indicates cooperativity of the network, and EC 50 , the stimuli dosage inducing half-maximal activity (representing threshold values for activation) for NF-κB activation obtained by experiment and simulation were 2.68 ± 6.92 and 1.48 × 10 − 6 ± 1.37 × 10 − 6 g/ml, 3.49 ± 0.168 and 6.07 × 10 − 7 ± 1.04 × 10 − 8 g/ml, respectively (Figures 1d and e). Furthermore, the simulation recapitulated the response in BCL10 and CARMA1 S578A mutants and wild-type cells in the presence of PKC and IKKβ inhibitions in our earlier study 13 Figure 1). We observed that the total sum of oscillation amplitudes was highly sensitive to perturbations of reactions involving IKKβ (66, 76, 77) and IκBα (142, 143; Figure 1j). The oscillation period was more sensitive to TAK1 (50), IKKβ (66) and IκBα (132, 140, 143; Figure 1k). Trends of parameter sensitivities for the sum of amplitudes and averaged period of NF-κB oscillations (Figures 1j and k) were generally consistent with those analyzed for individual peaks in the oscillations (Supplementary Figure 6). For the Hill coefficient and EC 50 values of NF-κB activity, highly sensitive parameters were mostly enriched in the upstream signaling reactions for TAK1 and IKKβ (48, 50, 65, 66, 76, 77; Figures 1l and m). These results show that reactions involved in positive feedback loops are strongly associated with oscillation amplitudes and periods. It is noteworthy that, although the switch-like induction of the first peak of nuclear NF-κB may significantly affect subsequent oscillation, the two dynamical behaviors are independently regulated by different components of the pathway.
Inhibition and enhancement of feedback loops The NF-κB network contains five feedback regulatory loops, two positive and three negative ( Figure 2). The positive feedback loops are fully contained in the upstream part of the network, whereas the negative feedback loops operate at the level or downstream of NF-κB (IκB's), or linking downstream and upstream network components (A20). It is known that the negative feedback loops mediated by transcriptional induction of IκBα and IκBε have distinct roles in controlling oscillation dynamics. 6,9,24 Consistent with previously reported results, we observed that changes in the parameters controlling NF-κB-induced IκBα expression (141) significantly altered oscillation amplitudes and periods ( Figure 2c). Changes in the parameter of NF-κB-induced IκBε expression (152) altered oscillation dynamics only at the late phase of the response (Figure 2d). We also determined that these negative feedback regulators had small effects on the Hill coefficient and EC 50 of NF-κB peak activity (Figures 2c and d).
On the other hand, the positive feedback loops operating upstream of NF-κB had more global effects. A stronger positive feedback from IKKβ to TAK1 (49) increased the amplitude of oscillations and the steepness of the switch responses ( Figure 2a). Changes in the parameter of the autoregulation of IKKβ (76) dramatically altered the oscillatory dynamics of NF-κB signals as well as the overall shape of the stimulus dose response ( Figure 2b). The negative feedback loop consisting of NF-κBinduced A20 expression (158) affected the Hill coefficient and EC 50 , as well as oscillation amplitudes and periods ( Figure 2e).
The above analysis predicted that the positive feedback loop from IKKβ to TAK1 would affect oscillation dynamics and the Hill coefficient (switch) but not the EC 50 (threshold) of the dose response of NF-κB activity ( Figure 2a). On the other hand, it was also predicted that changes in IκBα-related reactions would alter the oscillation patterns but not the switch-like response of NF-κB (Figures 1j-m). To test these predictions experimentally, we measured the time course and dose response of NF-κB activity using DT40 cells with either a CARMA1 S578 mutation to alanine (CARMA1 S578A ), which mimics the absence of positive feedback from IKKβ to TAK1, or IκBα overexpression (IκBα ox), which mimics strong negative IκBα feedback (  Table 3). As predicted, we found that in the CARMA1 S578A mutation, both oscillations and switch-like responses became indistinct (Figures 3a-e). The Hill coefficient of the CARMA1 S578A mutant (experiment; 0.286 ± 0.064, simulation; 0.530 ± 0.028) was lower compared with that in the wild-type cells, whereas the EC 50 did not change significantly (experiment; 3.74 × 10 − 7 ± 4.44 × 10 − 7 g/ml, simulation; 1.51 × 10 − 6 ± 2.55 × 10 − 7 g/ml; Figure 2a). In the IκBα overexpressing cells, the oscillations became indistinct, however, the steepness of the dose response did not change significantly (Figures 3f-j). The Hill coefficient and EC 50 values obtained by experiment (1.46 ± 0.204 and 2.57 × 10 − 7 ± 4.21 × 10 − 8 g/ml) and simulation (2.56 ± 0.124 and 6.69 × 10 − 7 ± 1.47 × 10 − 8 g/ml) for IκBα overexpression were not significantly different from those in wild-type cells. These observations indicate that oscillations and switch-like activation in the NF-κB system do not necessarily have to occur together, and can be separated.
Changes in protein abundance decouple oscillatory and switchlike responses To test the hypothesis that oscillations and switch-like responses can be separated, we next investigated the combinational effect of abundance of two proteins in the signaling pathway on the dynamics of NF-κB (Figure 4a; see numerical definitions of the dynamics in the 'Materials and Methods' section). Our results showed that most combinations of two proteins among BCL10, MALT1, CARMA1, TAK1, IKKβ and NF-κB result in both oscillations (red) and switch-like responses (blue) (Figure 4a, details are shown Switch and oscillation of NF-κB for B-cell activation K Inoue et al in Supplementary Figures 7 and 8). A notable exception to this trend was NF-κB. Our analysis predicts that a low abundance of NF-κB results in switch-like responses but not oscillations in some parameter regions, whereas high NF-κB abundance alone is likely to be sufficient to produce oscillations in other parameter regions. Specifically, the computational analysis predicted that cells with low MALT1 and high NF-κB abundances would show decoupling of oscillations and switch-like responses (Supplementary Figure 7). To confirm this experimentally, we determined NF-κB activity using MALT1 knockout (MALT1 − / − ) cells, with or without NF-κB Increased NF-κB abundance compensates for reduced upstream signaling activity Why do changes in NF-κB abundance cause oscillations in the absence of switch-like responses? To answer this question, we constructed a core model consisting of only core components, IKKβ, NF-κB and IκBα (Figure 5a, Supplementary Section 3), which are associated with critical parameters responsible for oscillations and switch-like responses (Figures 1j-m). Using this model, we performed equilibrium and quasi-equilibrium (also called nullcline) analysis (Supplementary Figure 9), which allowed us to qualitatively infer the type of dynamics expected near the equilibrium point.
The quasi-equilibrium analysis showed that oscillations are induced when the synthesis rate of the negative regulator IkB is transiently faster than its maximum degradation rate constant k 11 , which corresponds to IkBe in Figures 5b and c. Only an increase in NFkB* can overcome the degradation rate of IkB (yellow oscillation range, Figure 5b). An increase in totalNFkB causes an increase in NFkB* even if IKK* is low. Therefore, an increase in NF-κB abundance can induce adequate synthesis of IκB beyond its maximum degradation rate, even in spite of weak upstream signals, resulting in oscillations of NF-κB (black circles and green regions in Figures 5d-f). On the other hand, a decrease in totalIKK causes a major decrease in IKK* when the threshold is not exceeded (Supplementary Figure 9e), resulting in decreased or absent NFkB* oscillations (black circles and indigo regions in Figures 5d-f).
Roles of positive feedback in NF-κB oscillation dynamics and target gene expression We showed that the two positive feedback loops (IKKβ-TAK1 and IKKβ autoregulation) control the dynamics of NF-κB activity (Figures 2a, b and 3a-e, and Supplementary Figures 5m-o). Yet, protein abundance (Figure 4), equilibrium and quasi-equilibrium  Figure 5g). The time-course analysis showed that the first oscillation amplitude of NFkB* was much faster and higher in the presence of positive feedback and it was subsequently followed by prolonged oscillations (Figure 5h). IKK* also showed rapid velocity in the presence of feedback. This analysis demonstrated that the positive feedback loops enable prolonged oscillations of NF-κB activity by increasing the initial velocity of IKKβ and the amplitude of nuclear NF-κB. Prolonged oscillations in NF-κB activity are important for induction of target genes. 8,11,14 Therefore, we next elucidated the effect of the two responses of NF-κB dynamics on mRNA time course and dose response of CD83, a representative NF-κB target gene involved in B-cell maturation. 25 We tested four cell conditions; CARMA1 S578A mutant cells (no oscillation or switch response), IκBα overexpression cells (no oscillation but switch response), MALT1 mutant cells with NF-κB (RelA) overexpression (oscillation but no switch response) and wild-type cells (both oscillation and switch response; Figure 6). The CD83 gene expression level in the mutant cells was generally lower than in the wild-type cells, even though the NF-κB activity in MALT1 − / − /RelA ox cells was higher than in the wild-type cells (Figure 4d). This NF-κB overexpressing mutant, which is not associated with the switch response, also showed a slight increase in CD83 gene expression compared with the other mutants, implying that NF-κB oscillations themselves may have a small effect on NF-κB target gene expression. The wild-type cells also had a strong, prolonged and steep dose increase in CD83 gene expression. Interestingly, the Hill coefficient and EC 50 for CD83 gene expression in the wild-type cells was 8.19 ± 5.35 × 10 − 16 and 1.14 × 10 − 6 ± 2.93 × 10 − 8 g/ml (Supplementary Table 6), respectively, suggesting that cooperative

DISCUSSION
Positive feedback regulation has several functional roles in biological systems: to sustain molecular activity, to trigger oscillation or to assign a threshold, as seen, for example, in the MAPK and Cdc2/APC networks. 19,26 In the NF-κB network, the most important role of the positive feedback loops is to amplify the peak maxima of NF-κB nuclear translocation, which can induce the prolonged NF-κB activity necessary for target gene expression. Although two positive feedback loops are operative in the NF-κB signaling pathway, they do not provide sustained IKKβ activity to create prolonged NF-κB oscillation (Supplementary Figure 5d). Because the upstream IKK activity is transient, only increased NF-κB nuclear abundance can compensate for weak upstream IKK activity. Thus, the resulting high-amplitude NF-κB oscillations solely drive productive transcriptional responses important for activation of B cells. As a result of the current analysis, several key features were clarified as unique properties of the NF-κB network. First of all, switch-like activation and oscillatory behaviors of NF-κB are independently controllable through changes in kinetic parameters and abundance of proteins in the signaling network. Earlier studies indicated that NF-κB oscillation dynamics are mainly controlled by the IKK-IκB negative feedback loop. 6 However, our analysis showed that NF-κB oscillation is also controlled by stoichiometry and parameters associated with adaptor proteins such as CARMA1 and MALT1, molecules whose mutation and overexpression are responsible for the pathogenesis of some cases of B-cell lymphoma. 2,27,28 Furthermore, although trends of parameter sensitivities in the NF-κB network are quite similar for all oscillation peak amplitudes, the sensitivity values are increased by peak after the peak (Supplementary Figures 6a-d). This result suggests that oscillation amplitude, particularly in the late phase of oscillations, is more likely to be affected by small changes or noises, such as in external ligand concentrations and protein abundances in individual cells.
Our analysis provides rational evidence showing that NF-κB activity is under tight control of kinetics and abundance of proteins in the signaling pathway, and that dysregulation of this control can lead to abnormal NF-κB activation. In addition, switchlike expression of the CD83 gene suggests that the NF-κB network is a multifaceted, highly cooperative system for cellular commitment of B cells. Although we could not completely exclude the possibility of species differences, earlier studies strongly suggested that BCR signaling in chicken DT40 B cells can recapitulate, at least, the NF-κB activation pathway and its dynamics in primary mouse B cells; 13,29,30 therefore, our current study may shed light on the general mechanisms of gene expression regulated by NF-κB during B-cell differentiation.
In a theoretical aspect, because NF-κB shows damped oscillations, our mathematical analysis is different from the usual bifurcation analysis for attractors like sustained oscillations and switching phenomena between different periodic oscillations. So far, many mathematical studies on oscillatory dynamics in biological systems, such as cell cycle and circadian rhythms, have focused on sustained oscillations, 17,31 whose qualitative changes are tractable in terms of the conventional bifurcation theory. However, the mathematical formulation of qualitative changes in transient dynamics as presented in this study is highly significant as well, because such dynamics are responsible for the functionality of biological systems.

Generation of mutants and transfectants of DT40 cells
Chicken RelA and IκBα cDNA were generated by PCR. Each cDNA was cloned into the pApuro expression vector. 29 These constructs were transfected into MALT1-deficient 13 or wild-type DT40 cells by electroporation as described elsewhere. 29 Genomic clones of chicken RelA were obtained by PCR using oligonucleotides designed on the basis of NCBI database sequences (RelA; NM_205129) as primers and DT40 genomic DNA as a template. CARMA1 (wild type and its S578A mutant) and monomeric EGFP tagged wild-type RelA knock-in vectors were constructed as described previously 20 to generate a mini gene in the respective endogenous locus. Wild-type and various mutant DT40 cells were cultured in RPMI 1640 (GIBCO BRL, Carlsbad, CA, USA and Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum, 1% chicken serum, 50 μ mol/l 2-ME (Wako), 4 m mol/l L-glutamine and antibiotics.

Cell fractionation and western blots
For NF-κB activity, nuclear and cytoplasmic fractions were prepared as described previously. 13 In brief, the cells were lysed with lysis buffer containing 50 m mol/l Tris-HCl (pH 7.5), 0.5% Triton X-100, 137.5 m mol/l NaCl, 10% glycerol, 5 m mol/l EDTA and a proteinase inhibitor cocktail (Roche, Basel, Switzerland). The nuclei were separated by centrifugation and the pellets were resuspended in 1% NP-40 lysis buffer and used as the nuclear fraction. Western blot analysis was carried out after stimulation with anti-IgM (M4) as described previously. 13,20,29 The samples were subjected to western blot, and RelA activity was quantified from the intensities of protein bands using a Multi Gauge version2.2 (Fujifilm, Tokyo, Japan) densitometer. 13 For protein detection, the ECL Plex fluorescent western blotting system and ImageQuant LAS 4000 (GE Healthcare, Little Chalfont, UK) were used.

Mathematical model
A comprehensive diagram of the BCR-activated NF-κB signaling network is shown in Supplementary Figure 1. The mathematical model was built on the basis of 47 ordinary differential equations with 196 parameters. Details about the rate laws and parameter values can be found in Supplementary Table 1. Steady state values in resting cells were obtained to simulate the basal level (initial values) of the cell before BCR stimulation. Then, computer simulation was performed using the function for the input signal shown in Supplementary Figure 2a. Simulation was carried out with the function CVODE (http://computation.llnl.gov/casc/sundials/main.html) in C language. Parameter values were searched by genetic algorithm with Unimodal Normal Distribution Crossover 32 and Minimal Generation Gap. 33 In genetic algorithm, we used the following cosine similarity as fitness function.
where CS is cosine similarity, n is the number of molecules (CARMA1, TAK1, IKK or NF-κB), w i is a weight parameter, x ! simulation;i is vector for values of the molecule i at each time point in the simulation, x ! experiment;i is vector for activities of the molecule i at each experimental time point, is the inner product of x ! simulation;i and x ! experiment;i , and U j j denotes the magnitude of a vector. A low cosine similarity means that dynamic behaviors of the molecules in simulation are similar to the ones in experiment. We used the parameter set with the lowest cosine similarity in our parameter search.

Definition of amplitudes and periods of oscillation
Quantitative measures for amplitudes and periods of oscillation of NF-κB (Figures 1f and g) were defined as follows: the first peak amplitude is defined as the difference between the height of the first peak and the basal value before cell stimulation. The second amplitude is defined as the difference between the height of the second peak and the value of the bottom after the first peak. The third and fourth amplitudes were defined in the same way as described above. Periods were defined as the time lengths between the peaks.

Hill equation
The Hill equation was defined as follows: where 'bottom' is minimum activity, 'top' is maximum activity, EC50 is halfmaximum effective dose, h is Hill coefficient, y is NF-κB activity and x is stimulus dose. The parameters bottom, top, h and EC50, were optimized using the function NonLinearModel.fit of MATLAB R2014a (MathWorks, Natick, MA, USA).

Sensitivity analysis
The single parameter sensitivity of each reaction is defined by where v i is an ith reaction, v is reaction vector v = (v 1 , v 2 , …), and q(v) is a target function, e.g., amplitude, period, Hill coefficient, EC 50 . The total sum of amplitudes of all peaks, and the average of periods were also used for the analysis. The sensitivity was calculated with 1% increases in the reaction rates.
Numerical criteria for oscillation and switch-like response Numerical criteria of oscillations and switch-like responses were defined as follows: oscillations have more than two peaks, which have relative amplitude values of 10 − 2 to the maximum activity. Criteria of switch-like dose responses were defined as one satisfying the following three conditions; Hill coefficient (41), EC 50 ( o10 − 4 ), and dose response ratio of more than 2-fold between maximum and minimum activity values.