Email updates

Keep up to date with the latest news and content from BMC Systems Biology and BioMed Central.

Open Access Research article

Interactions among oscillatory pathways in NF-kappa B signaling

Yunjiao Wang14*, Pawel Paszek2, Caroline A Horton2, Douglas B Kell3, Michael RH White2, David S Broomhead4 and Mark R Muldoon4*

Author Affiliations

1 Mathematical Biosciences Institute, The Ohio State University, Jennings Hall, Columbus, Ohio 43210, USA

2 Centre for Cell Imaging, School of Biological Sciences, Bioscience Research Building, Crown Street, Liverpool, L69 7ZB, UK

3 School of Chemistry and The Manchester Interdisciplinary Biocentre, University of Manchester, 131 Princess Street, Manchester M1 7DN, UK

4 School of Mathematics, Alan Turing Building, University of Manchester, Manchester M13 9PL, UK

For all author emails, please log on.

BMC Systems Biology 2011, 5:23  doi:10.1186/1752-0509-5-23

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1752-0509/5/23


Received:26 April 2010
Accepted:3 February 2011
Published:3 February 2011

© 2011 Wang et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

Sustained stimulation with tumour necrosis factor alpha (TNF-alpha) induces substantial oscillations—observed at both the single cell and population levels—in the nuclear factor kappa B (NF-kappa B) system. Although the mechanism has not yet been elucidated fully, a core system has been identified consisting of a negative feedback loop involving NF-kappa B (RelA:p50 hetero-dimer) and its inhibitor I-kappa B-alpha. Many authors have suggested that this core oscillator should couple to other oscillatory pathways.

Results

First we analyse single-cell data from experiments in which the NF-kappa B system is forced by short trains of strong pulses of TNF-alpha. Power spectra of the ratio of nuclear-to-cytoplasmic concentration of NF-kappa B suggest that the cells' responses are entrained by the pulsing frequency. Using a recent model of the NF-kappa B system due to Caroline Horton, we carried out extensive numerical simulations to analyze the response frequencies induced by trains of pulses of TNF-alpha stimulation having a wide range of frequencies and amplitudes. These studies suggest that for sufficiently weak stimulation, various nonlinear resonances should be observable. To explore further the possibility of probing alternative feedback mechanisms, we also coupled the model to sinusoidal signals with a wide range of strengths and frequencies. Our results show that, at least in simulation, frequencies other than those of the forcing and the main NF-kappa B oscillator can be excited via sub- and superharmonic resonance, producing quasiperiodic and even chaotic dynamics.

Conclusions

Our numerical results suggest that the entrainment phenomena observed in pulse-stimulated experiments is a consequence of the high intensity of the stimulation. Computational studies based on current models suggest that resonant interactions between periodic pulsatile forcing and the system's natural frequencies may become evident for sufficiently weak stimulation. Further simulations suggest that the nonlinearities of the NF-kappa B feedback oscillator mean that even sinusoidally modulated forcing can induce a rich variety of nonlinear interactions.

Background

Nuclear factor kappa B (NF-κB) transcription factors are critical to the control of response to cellular stress and are also involved in the regulation of cell-cycle/growth, survival, apoptosis, inflammation and immunity [1-5]. They are dimeric molecules, composed of either homo-or hetero-dimers, with the most common form being the RelA:p50 hetero-dimer.

In resting cells, NF-κB dimers are sequestered in the cytoplasm by members of a family of molecules, the inhibitors of κB or IκBs. When the cells are stimulated with tumour necrosis factor alpha (TNFα) certain kinases, the IKKs, are activated and phosphorylate IκBs, marking them for degradation. The degraded IκBs then release NF-κB dimers which are free to translocate into the nucleus, where they bind to specific sequences in the promoter or enhancer regions of target genes, including those for IκBα and the zinc finger protein A20 [6,7]. Newly synthesized IκBα migrates to the nucleus, binds to NF-κB dimers and removes them from the nucleus, while A20 protein stays in cytoplasm and represses the activity of TNFα receptors [8]. Hence the NF-κB system includes at least two negative feedback loops, one involving cytoplasmic sequestration mediated by IκBα and another involving A20.

Coupled negative feedback loops are known to have the potential to support oscillations [9,10]. Biological evidence suggesting NF-κB oscillations was reported by Hoffmann et al. [11], who did population-level studies of IκBα-/-embryonic fibroblasts using electro-mobility shift assays, and also observed by Nelson et al. [12] in the cellular concentration and nuclear: cytoplasmic localisation of IκBα and RelA-fluorescent fusion proteins in single SK-N-AS cells. Single cell time-lapse imaging data showed persistent cycling of NF-κB localisation between the cytoplasm and nucleus of SK-N-AS cells in response to continual TNFα stimulation [12,13].

The expression of genes regulated by NF-κB is tightly coordinated with the activities of many other signalling and transcription-factor pathways [14-18] including the p53 signalling pathway. Though the canonical NF-κB signalling pathway has been studied extensively, the existence and mechanisms of the interactions between the NF-κB pathway and other signalling pathways are unclear.

Ashall et al. [19] reported that different patterns of pulsatile stimulation of the NF-κB system lead to different patterns of NF-κB dependent gene expression, supporting the view that the frequency of the oscillations may have a functional role. It is therefore important to understand how the system's response frequencies are influenced by interactions with other oscillatory pathways. If both the NF-κB network and its couplings to other other oscillatory pathways were purely linear, then it would be straightforward to use the machinery of transfer functions to characterise their interactions. In particular, one would expect the power spectrum of a periodically stimulated system to have its power concentrated at the forcing signal's frequency and its harmonics. But the NF-κB network is a highly nonlinear system of coupled chemical reactions and hence, as we will demonstrate in modelling studies below, its response to periodic stimulation can include complex interactions between the intrinsic negative feedback oscillator and the stimulus. In what follows we work with a deterministic model first described in [19] and examine the power spectral densities of the time courses of NF-κB localisation when the model is subjected to two types of periodic signals: trains of rectangular pulses and sinusoidal signals. Both sorts of stimulus are intended as a proxies for the influence of other oscillatory pathways on the core NF-κB feedback loop and we find that such interactions can produce a rich variety of nonlinear dynamical behaviour. Our results on sinusoidal forcing are in broad agreement with those of Fonslet et al. [20], who applied a related family of sinusoidal stimulation protocols to a reduced model of the core NF-κB oscillator originally developed in [21].

Results and Discussion

Experimental phenomena and time series analysis

We begin by analyzing the time courses of the localisation of fluorescently tagged NF-κB in single cells subjected to TNFα stimulation (see [22] for details). The panels in the top row of Figure 1 show the ratio of nuclear-to-cytoplasmic fluorescence for three patterns of stimulation. In all cases the cells received three strong (10 ng/ml) pulses of TNFα, each of five minute duration. These pulses were separated by various inter-pulse intervals—55 minutes for panel (a),95 minutes for panel (b) and 195 minutes for panel (c)—resulting in net stimulus periods of 60, 100 and 200 minutes. The first of these is rather shorter than the observed period of the oscillations induced by constant TNFα stimulation, which is close to 100 minutes.

thumbnailFigure 1. Power spectra of single-cell data. Panels (a), (b) and (c) show time courses of the ratio of nuclear-to-cytoplasmic fluorescence recorded from single cells while panels (d), (e) and (f) show the corresponding power spectral densities, normalized to have unit area beneath the curve. Note that the frequencies in panels (d)-(f) are measured in units of the stimulus frequency and so the horizontal scale varies from panel to panel.

The bottom row of Figure 1 shows the corresponding power spectral densities and makes the point that the intrinsic oscillation—the one induced by constant stimulation—does not appear to have been excited by these pulsed stimulation protocols: in all cases the peaks in the power spectral density appear at multiples of the stimulus frequency. In mathematical terms, the cells respond as though they have been perturbed away from a stable resting state and do not offer evidence of any more complex internal dynamics. Nonetheless, these data provided crucial input to the development of the models in [19]. In particular, note that the amplitude of the response to the second and third pulses of stimulation varies as a function of the inter-pulse period, being somewhat reduced in panels (a) and (b), but not in (c): this constrains the rate at which the model should relax toward equilibrium.

One possibility is that the pulsing signal is too strong and suppresses the response at the natural frequency. In the following section we report in numerical experiments, based on the deterministic model in [19], through which we investigate this possibility by simulating the system's response to a periodic train of square pulses with strength varying from 0 to 10 ng/ml and inter-pulse intervals ranging from 20 to 295 minutes.

Model Introduction and Bifurcation Analysis

In this section, we introduce a deterministic mathematical model developed in [19]. It includes two coupled negative feedback loops—one for NF-κB-IκBα interactions and another that models the dynamics of A20—as is illustrated in Figure 2. Here A20 acts by modifying the activity of IκB kinase (IKK), which transduces the TNFα signal by phosphorylating IκBα, thus initiating the production of free NF-κB. IKK is assumed to exist in one of three states: neutral IKK (denoted IKKn, a state ready for activation), active IKK (the state capable of phosphorylating IκBα) and inactive IKK (IKKi, a state unable to phosphorylate IκBα and also incapable of being activated by the TNFα signal). TNFα activates IKK by transforming neutral IKK into active IKK, which then acts to liberate NF-κB. Active IKK is assumed to convert to its inactive form spontaneously, with linear kinetics. In the absence of A20, this inactive IKK would then transform back into neutral IKK, also with linear kinetics. A20 down-regulates NF-κB activity indirectly, by retarding the transformation of inactive IKK to the neutral form: Equations (1) provide a simple model for these dynamics:

thumbnailFigure 2. The deterministic, two-loop model of [19]. Proteins are shown as coloured ovals, while mRNA's are shown round-cornered rectangles. In addition to the direct negative feedback through which nuclear NF-κB leads to the production of IκBα, and thus to its own inactivation, the model also includes an indirect negative feedback mediated by A20's influence on the signal transduction machinery.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M1">View MathML</a>

(1)

Here we have suppressed the time dependence of the concentrations of the various forms of IKK and of the A20 protein. The k* are fixed parameters and TR is a parameter indicating the presence or absence TNFα stimulation: TR = 1 when the system is being stimulated with 10 ng/ml TNFα and TR = 0 otherwise. This is a simplification: TNFα does not act directly on IKK, but rather produces its effect through a cascade of chemical reactions that begins when TNFα binds its receptor at the cell surface.

Following Ashall et al., we treat only the final stages of this chain. Thus when we model strong stimulation by setting TR = 1, we abstract away a great deal of transduction machinery. In the studies that follow we will want to explore the consequences of weaker stimulation (TNFα concentrations of 10 ng/ml are orders of magnitude higher than physiological levels: Matalka et al. [23] report measurements from various tissues in healthy mice and found concentrations in the range 1-5 pg/ml while Prabha et al. [24] studied the plasma of healthy human subjects and found TNFα concentrations on the order of 100 pg/ml.) and it is natural to generalize the role of TR, allowing it to vary across the interval 0 ≤ TR ≤ 1. In light of the modelling assumptions discussed above, one should not imagine that TR depends linearly on the concentration of TNFα, but only that TR increases monotonically with dose.

The model of Ashall et al. captures two important dynamical features observed in single cell data [12,19]: (i) continuous stimulation with high doses of TNFα leads to sustained oscillations with a period of around 100 minutes and (ii) pulsatile stimulation with the same high concentration of TNFα leads—as illustrated in Figure 1—to entrainment of the response by the pulsing signal. Given that continuous, strong stimulation induces sustained oscillations, one is prompted to ask how the existence of these oscillations depends on the strength of stimulation. We addressed this question by doing a bifurcation analysis using TR as the parameter: Figure 3 is the resulting diagram. The model has a Hopf bifurcation (HB) at TR = T* ≈ 0.366, which means that for TR > T* the ratio of nuclear-to-cytoplasmic NF-κB concentration exhibits sustained oscillations, but for TR < T* only damped oscillations. A forthcoming paper, Wang et al. [25], will present a comprehensive survey of the bifurcation structure of this model, but here we will concentrate on spectral analysis of the responses.

thumbnailFigure 3. Qualitative behaviour as a function of TR. At left, a bifurcation diagram showing how the oscillatory response depends on TR. The solid black curve in the region 0 ≤ TR T* represents a branch of stable equilibria-steady, non-oscillating solutions-while the dashed black curve for TR >T* indicates a branch of unstable steady states. The red curves show the limits-the peak and trough values-of the stable oscillatory responses that exists for these values of TR. The panel at right shows period of the oscillation as a function of TR.

Numerical experiments

All the numerical experiments reported here begin from the same initial condition, whose preparation is described in the Methods section. We simulated the consequences of two TNFα stimulation protocols: pulse-like stimulation similar to that used in Nelson's experiments and sinusoidally modulated stimulation.

Pulsed Stimulation

From the experiments and simulations in [19], we know that stimulation with a finite train of three strong pulses tends to entrain the cell's response. In this section, we first show that when subjected to a long periodic train of strong pulses, the model's response frequencies are also entrained, but that when the model is driven by sufficiently weak periodic pulse trains, various resonance phenomena connected to the underlying natural oscillations become observable.

Figure 4 illustrates the continuously-pulsed analogue of the experiments from Figure 1. It shows the spectral content of the model's response to a long periodic pulse train in which five minute periods of strong (TR = 1) stimulation alternate with periods during which TR = 0. The panel at left shows a typical response to this sort of strong, pulsed, periodic forcing: after a brief transient the response is periodic with the same period as the forcing pulse train. The right panel, which is a heat map showing the power spectral density of the response as a function of pulsing frequency, shows that similar entrainment occurs over a wide range of frequencies: the bright lines—which correspond to peaks in the power spectral density—have integer slopes, indicating that the power in the response is concentrated at harmonics of the forcing frequency.

thumbnailFigure 4. Periodic stimulation with strong pulses. The panel at left shows the model's response to sustained periodic stimulation in which strong (TR = 1) pulses of 5 minute duration alternate with stimulus-free intervals of 55 minute duration. The panel at right shows the power spectral density (PSD) of the response as a function of pulsing frequency: a vertical strip cut from this panel provides a colour-coded version of the PSD of the models's response to pulsed forcing at the corresponding frequency.

However, when the stimulating pulses are weaker the nonlinearity of the system leads to more complex patterns of resonance. In particular, both subharmonic resonance—when periodic forcing excites a response at a rationally-related lower frequency—and superharmonic resonance—in which a pure sinusoid excites responses containing higher harmonics—are possible. Note that these terms are defined with reference to the forcing frequency, a convention used in, for example, [26]. Superharmonic resonance is harder to identify when, as with the rectangular pulses used here, the periodic forcing already has power at higher harmonics, but subharmonic resonance occurs when 0.01 < TR < 0.2, as is evident in the power spectral densities summarized in Figure 5. In addition to the bright lines with integer slopes, lines with slopes <a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M2">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M3">View MathML</a> also appear. The first of these provides evidence that subharmonic resonance occurs over a wide range of pulsing frequencies.

thumbnailFigure 5. Response to periodic trains of weak pulses. A heat map summarising the power spectral density of NF-κB localisation when the model is subjected to pulsed forcing with TR = 0.1. For forcing frequencies in the range 1.5-3.0 × 10-4. Hz. the lowest-lying line segment, which has slope ≈ 1/2, provides evidence of subharmonic resonance.

When TR ≥ 0.2 the power spectral density of the response shows power only at integer multiples of the pulsing frequency, indicating that the response is fully entrained by the forcing. Given that the pulse-strength used in the experiments corresponds to TR = 1, these modelling results are in qualitative agreement with the experimental data in Figure 1: the stimulation was so strong that we should have expected it to have entrained the responses completely.

Sinusoidally modulated stimulation

In this section we study the response of Horton's model to sinusoidally modulated stimulation of the form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M4">View MathML</a>

(2)

Here 0 ≤ ε ≤ 1/(1 + η) is the time average of the stimulus strength while 0 ≤ η ≤ 1 and ν > 0 are, respectively, the relative amplitude and the frequency of the sinusoidal modulation. With this parameterization TR (t) has period τ = 1/ν and range

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M5">View MathML</a>

The roles of the parameters are illustrated in Figure 6.

thumbnailFigure 6. Sinusoidal forcing function. The function TR (t) for the parameter values ν = 1, η = 0.5 and ε = 0.6.

The case with η = 1 is the straightforward substitution of pulse trains with sinusoidal waves having peak-to-trough amplitude 2ε. Our motivations for this formulation are twofold: firstly, provided η < 1 we have TR(t) > 0 at all times and so expect the system more readily to exhibit oscillations similar to those induced by constant stimulation. Secondly, periodic forcing of the form (2) converts our model into a sinusoidally forced nonlinear oscillator, a class of systems that has been studied very extensively (see, for example, Pikovsky et al. [27], Wiggins [26] or Nayfeh and Mook [28]).

When η = 0 the forcing (2) reduces to a constant stimulus with TR (t) = ε and the bifurcation analysis illustrated in Figure 3 leads us to expect stable oscillations when ε > 0.366 and a stable steady-state otherwise. But for η > 0 the system's response depends delicately on the relationship between the forcing frequency ν and the natural frequency ν0. When the forcing frequency ν and natural frequency ν0 are close (the requisite degree of closeness depends on η) the response will become mode locked, or synchronized with the forcing: it will then have the same frequency ν as the forcing and its power spectral density will be concentrated around the harmonics of ν. When the difference (ν -ν0) is larger—when it lies just outside of a critical interval around zero whose size depends on η —the response becomes quasiperiodic and its power spectrum has features at frequencies f of the form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M6">View MathML</a>

(3)

where p and q are integers.

If and 0 are close, the system's nonlinearity will permit what is called synchronization of order p/q: there will be periodic responses in which the intrinsic oscillator goes through p cycles for every q periods of the forcing, so that the response has frequency

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M7">View MathML</a>

(4)

This idea allows us to give precise definitions for the terms sub-and superharmonic resonance used above: the former corresponds to resonances where p = 1 and q > 1 in (4), while the latter corresponds to q = 1 and p > 1. Finally, when both η and the difference (ν -ν0) are large, the system's response can become chaotic, so that the features in its power spectral density bear no simple relationship to the frequencies ν and ν0.

Figure 7 which is a heat map of the power spectral density of the response generated by relatively strong sinusoidal forcing of the from (2), illustrates many of the behaviours described above. Consider first the vertical strip with ν/ν0 ≈ 1. In this region the sinusoidal modulation entrains the NF-κB system essentially completely and so the heat map resembles the corresponding region in the left panel Figure 4: the power in the response is concentrated along lines corresponding to harmonics of the forcing frequency. Away from the region ν/ν0 ≈ 1 the power spectrum of the response is considerably more complex: the strong horizontal bands in Figure 7 which occur at integer multiples of ν0, show that there is substantial power at the NF-κB system's natural frequency and its harmonics. Additionally, the nonlinearity of the system means that the response has power at frequencies given by

thumbnailFigure 7. Nonlinear resonances to sinusoidal forcing. A heat map showing the power spectral density of the response of Horton's model to forcing of the form (2) with ε = η = 0.5 and frequencies in the range 0 ≤ ν ≤ 4.5 × ν0.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M8">View MathML</a>

which gives rise to the network of lines with slope of ±1. Finally, Figure 7 also exhibits sub-harmonic resonance: vertical strips for which ν 0 (with q a whole number) show power concentrated at the forcing frequency, but also at frequencies f = ν/p, where p is a whole number. This is perhaps clearest in the strip ν ≈ 2ν0, where the strongest spectral feature lies along the line f = ν/2 ≈ ν0.

The two panels of Figure 8 are analogues of Figure 7 but with especially strong (η = 1, left panel) or weak (η = 0.1, right panel) modulation. The qualitative features are much the same, though it is interesting to note that in the limit of very strong modulation—when η = 1 and so TR (t) vanishes once per forcing period—the system is very strongly entrained by the forcing and does not show much power at its natural frequency ν0 or its harmonics until the forcing frequency ν > 1.5 ν0. By contrast, when the modulation is weak evidence of modal interaction is also weak, with very narrow mode-locking regions near forcing frequencies of the form ν 0.

thumbnailFigure 8. Weak and strong modulations. Heat maps showing the power spectral density of the responses to forcing of the form (2) with ε = 0.5 and η = 1 (left) or η = 0.1 (right).

As Figures 7 and 8 illustrate, the susceptibility of our model NF-κB system to resonance with sinusoidal modulations depends strongly on the amplitude of the modulation: Figure 9 provides a quantitative survey of this phenomenon. The shaded regions are examples of Arnol'd tongues: their precise shapes can be calculated with the methods outlined in the Appendix. When the modulation frequency and amplitude (ν, η) lie inside these tongues, the response of the system will be periodic, with a period τ that is an integer multiple of the modulation period 1/ν and that also lies close to an integer multiple of the natural Period 1/ν0:

thumbnailFigure 9. Arnol'd tongues. The qualitative behaviour of Horton's NF- κB system when subjected to forcing of the from 2 with ε = 0.5: when the modulation frequency ν and amplitude η lie inside the shaded regions the response is periodic with a period that is related, by a relation like (5), to both the modulation's period and the system's natural period.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M9">View MathML</a>

(5)

where p and q are whole numbers. Generally speaking these periodic responses—which are a nonlinear generalization of the familiar phenomenon of resonance in linear systems—are easiest to excite when the numbers p and q are small: the tongues in Figure 9 are labelled by ratios q:p where p and q are as in (5). Although Figure 9 illustrates this story for the specific family of modulations with ε = 0.5, the qualitative picture is essentially the same for all values in the range of 0.366 < ε < 0.5: in all cases there will be narrow tongues of parameter combinations (ν, η) for which the response is periodic with a period given by a resonance relation like (5). For modulations whose parameters lie outside the tongues the qualitative behaviour of the response will be more complex, having a power spectrum qualitatively similar to those in Figure 7, with features at harmonics of the modulation frequency, at harmonics of the natural frequency ν0 (which varies with ε) and at near-resonant combinations (3) of the two.

When TR < 0.366, the unforced system is a damped oscillator and so, when subjected to periodic forcing, may exhibit resonant phenomena. By using the same numerical simulation and power spectral density analyses as we did for the case ε = 0.5, we find that the responses of the pathway to forcing with with TR < 0.366 can be divided into two groups: those in which harmonic, sub-and superharmonic resonances can be observed (in the region 0 < ε ≤ 0.01), and those for which only harmonic and subharmonic resonance can be observed.

We have focussed here on periodic and quasiperiodic oscillatory interactions amenable to power-spectral analyses of the sort illustrated in Figures 4, 5, 6, 7, 8 and, for single-cell recordings, in the lower panels of Figure 1. But there is every reason to expect a much richer range of dynamical behaviour: Fonslet et al. [20], who applied forcing of the form (2) with η = 1 to a simplified NF-κB model, saw evidence of a period-doubling cascade as well as chaos and strange attractors. Complete analysis of these more complex dynamical regimes requires tools beyond the power spectrum, and so we defer their exploration to a future paper.

Conclusions

Given that intrinsically nonlinear chemical kinetics underpin cell-signalling networks, one shouldn't expect these systems to be linear and any analysis of experimental time series, whether in the time or frequency domain, must take this intrinsic nonlinearity into account. Our modelling studies suggest that coupling even the simplest, sinusoidal signal into the NF-κB network can give rise to a host nonlinear phenomena, including harmonic, subharmonic and superharmonic resonances as well as quasiperiodic and even chaotic behaviour.

The simulation studies reported here used an externally-imposed oscillatory forcing as a proxy for interactions between the core NF-κB feedback loop and other oscillatory networks. Our results suggest that interactions between the NF-κB oscillator and other oscillatory pathways can give rise to extremely rich temporal signalling programs and so, perhaps, to many distinct patterns of expression for target genes.

Methods

Differential Equations

The model studied in this paper is specified by the following ODEs. External forcing via TNFα stimulation is represented by the function Tr (t), whose value ranges between 0 and 1, indicating the strength of stimulation as a fraction of the maximum possible.

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M10">View MathML</a>

Initial conditions

All the numerical experiments reported above began with the initial concentrations

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M28">View MathML</a>

while all others start at zero. We then integrated the ODEs through 2000 minutes of simulated time and used the result as an initial condition for the various TNFα-stimulation studies.

Numerics

The bifurcation diagram in Figure 3 and the plot of the Arnol'd tongues in Figure 9 were prepared using Bard Ermentrout's XPPAUT [29], a front end for the powerful numerical bifurcation package AUTO [30]. All our other numerical work used MATLAB© to integrate the ODEs, compute the power spectra and plot the figures.

Authors' contributions

YW did the numerical experiments reported here, basing her work on code originally written by CAH, who also developed the underlying model. PP and CAH provided biologically-informed advice on modelling, especially on how best to couple the model to an external signal. MRHW, DBK and DSB provided critical readings of the manuscript, as well scientific guidance throughout the project. MRM worked closely with YW on the computational side and, with YW, wrote the first draft and prepared the figures. All authors reviewed and approved the final draft.

Appendix: circle maps

Here we discuss a standard mathematical tool, the circle map, used to study the response of a nonlinear oscillator subjected to periodic forcing. The idea is to consider sufficiently weak forcing that the response is close to that of the unforced system, and then define a phase angle ϕ that is close to the phase of the corresponding unforced oscillator.

If ω0 is the frequency of the unperturbed oscillator and ω is the forcing frequency then, after one period of the forcing, ϕ will have advanced by

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M11">View MathML</a>

where τ = 2π/ω is the period of the forcing, ε is a measure of its amplitude and G(ϕ) is a phase-dependent function that characterizes the response.

The simplest circle map, studied here by way of illustration, is given by a function F : [0, 2π) → [0, 2π) of the form of

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M12">View MathML</a>

(6)

where α = (ω0/ω), and 0 < ε ≪ 2π. The dynamics of a circle map can be characterized with a single parameter called the rotation number. For a given initial point ϕ0, the rotation number is defined as the long-term average phase shift per one iteration (that is, per one period of the forcing):

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M13">View MathML</a>

(7)

where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M14">View MathML</a>

One can show that the limit in the definition of the rotation number (7) exists and does not depend on the initial point ϕ0 [31], p.102. The qualitative dynamics of the forced system are thus of two types: motions with rational rotation numbers and those with irrational ones. Further, the rotation number is retional-say <a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M15">View MathML</a>, with p ∈ ℕ and q ∈ ℤ+ —if and only if there exists some ϕ0 ∈ [0, 2π) such that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M16">View MathML</a>

(8)

That is, rational rotation numbers correspond to periodic dynamics: the orbit of ϕ0 is called a (p, q) cycle. On the other hand, irrational rotation numbers correspond to quasiperiodic dynamics [27,31,32].

According to a theorem of Denjoy [33], if the rotation number is irrational there is a continuous, invertible change of coordinates h: [0, 2π) → [0, 2π), say, h(ϕ) = θ, such that

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M17">View MathML</a>

Where <a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M18">View MathML</a>. This implies that forced oscillators whose corresponding circle maps have irrational rotation number never repeat periodically.

Results of Arnold [31] show that the qualitative behaviour of the circle map (6) is stable against arbitrary small perturbations if and only if the rotation number is rational. We can thus expect that if F has a rational rotation number <a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M19">View MathML</a>, there exists a region of parameter values (α, ε) such that the all the forced systems whose parameters lie in this region share the same rational rotation number: that is, they all have the same sorts of periodic response. Such regions of parameter values are called Arnol'd tongues.

A similar story holds for the sinusoidally forced NF-κB system: when the time-average of the stimulation is sufficiently strong that the system supports periodic oscillations with frequency ω0, and when the amplitude of the forcing is sufficiently weak that it is sensible to measure a phase ϕ with respect to that of the intrinsic oscillator, then one can construct numerically a circle map of the form

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M20">View MathML</a>

(9)

where T = 2π/ω is the period of the forcing. Although the function G' in (9) is not as simple as the the G (ϕ) in (6), the map F' still has (p, q) periodic responses and a corresponding system of Arnol'd tongues: Figure 9 shows examples.

We'll conclude this Appendix by explaining why, at least for the simple circle map (6), the Arnol'd tongues are wedge-shaped. Consider the case where F has rational rotation number <a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M21">View MathML</a>. Then, as mentioned above, there exists a periodic point ϕ0 satisfying (8). If we expand Fq (ϕ0) in powers of ε we find

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M22">View MathML</a>

Consider the case where α is very close to a rational number: <a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M23">View MathML</a> with |β| ≪ 1/q. Then

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M24">View MathML</a>

Where

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M25">View MathML</a>

(10)

It is not hard to show that GP/q is bounded and continuous when regarded as a function of ϕ0 and so attains its maximum and minimum values on [0, 2π]. Thus, for each fixed ε, there is an interval of β on which GP/q = 0 for some ϕ ∈ [0, 2π]. Now consider the way in which the end points of this interval depend on ε. The case q = 1 and p = 0 is especially clear: Eqn. (10) becomes

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M26">View MathML</a>

which has some solutions ϕ ∈ [0, 2π] provided that |2πβ| ≤ ε or |β| ≤ ε/2π.

The values q = 1 and p = 1 correspond to the case where the periodic response has the same frequency as the forced system and thus undergoes one complete phase-rotation per period of the forcing. The analysis sketched above shows that this particular form of mode-locking will pertain for (ε, α) in a wedge-shaped region whose width increases linearly with forcing amplitude ε. For more general mode-locking resonances, where q ≥ 2, Arnol'd showed in [31] that the width of the tongue satisfies:

<a onClick="popup('http://www.biomedcentral.com/1752-0509/5/23/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1752-0509/5/23/mathml/M27">View MathML</a>

Acknowledgements

CAH and YW were supported by BBSRC grant BBD0088081. PP and MRHW were supported by MRC Grant G0500346. PP, MRM, DSB and MRHW were supported by BBSRC grants BBF005381 & BBF00561X1.

References

  1. Hayden M, Ghosh S: Signaling to NF-κB. [http://dx.doi.org/10.1101/gad.1228704] webcite

    Genes & Development 2004, 18:2195-2224. OpenURL

  2. Bonizzi G, Karin M: The two NF-κB activation pathways and their role in innate and adaptive immunity.

    Trends in Immunology 2004, 25(6):280-288. PubMed Abstract | Publisher Full Text OpenURL

  3. Pasparakis M, Luedde T, Schmidt-Supprian M: Dissection of the NF-κB signalling cascade in transgenic and knockout mice.

    Cell Death & Differentiation 2006, 13:861-872. OpenURL

  4. Pahl HL: Activators and target genes of Rel/NF-κB transcription factors. [http://www.nature.com/onc/journal/v18/n49/abs/1203239a.html] webcite

    Oncogene 1999, 18(49):6853-6866. PubMed Abstract | Publisher Full Text OpenURL

  5. Ghosh S, May M, Kopp E: NF-κB and rel proteins: Evolutionarily conserved mediators of immune responses.

    Annual Review of Immunology 1998, 16:225-260. PubMed Abstract | Publisher Full Text OpenURL

  6. Scott ML, Fujita T, Liou HC, Nolan GP, Baltimore D: The p65 subunit of NF-kappa B regulates I kappa B by two distinct mechanisms. [http://genesdev.cshlp.org/content/7/7a/1266.abstract] webcite

    Genes & Development 1993, 7(7a):1266-1276. OpenURL

  7. Song HY, Rothe M, Goeddel DV: The tumor necrosis factor-inducible zinc finger protein A20 interacts with TRAF1/TRAF2 and inhibits NF-kappaB activation.

    Proceedings of the National Academy of Sciences of the United States of America 1996, 93(13):6721-6725. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. Wertz IE, O'Rourke KM, Zhou H, Eby M, Aravind L, Seshagiri S, Wu P, Wiesmann C, Baker R, Boone DL, Ma A, Koonin EV, Dixit VM: De-ubiquitination and ubiquitin ligase domains of A20 downregulate NF-κB signalling.

    Nature 2004, 430:694-699. PubMed Abstract | Publisher Full Text OpenURL

  9. Pigolotti S, Krishna S, Jensen MH: Oscillation patterns in negative feedback loops.

    Proceedings of the National Academy of Sciences 2007, 104(16):6533-6537. Publisher Full Text OpenURL

  10. Thomas R: Logical analysis of systems comprising feedback loops.

    Journal of Theoretical Biology 1978, 73(4):631-656. PubMed Abstract | Publisher Full Text OpenURL

  11. Hoffmann A, Levchenko A, Scott M, Baltimore D: The I kappa B-NF-kappa B signaling module: Temporal control and selective gene activation.

    Science 2002, 298(5596):1241-1245. PubMed Abstract | Publisher Full Text OpenURL

  12. Nelson DE, Ihekwaba AEC, Elliott M, Johnson JR, Gibney CA, Foreman BE, Nelson G, See V, Horton CA, Spiller DG, Edwards SW, McDowell HP, Unitt JF, Sullivan E, Grimley R, Benson N, Broomhead DS, Kell DB, White MRH: Oscillations in NF-kappa B signaling control the dynamics of gene expression.

    Science 2004, 306(5696):704-708. PubMed Abstract | Publisher Full Text OpenURL

  13. Horton CA: Computational Modelling of Cell Signalling Pathway Dynamics. PhD thesis. School of Biological Sciences, University of Liverpool; 2006. OpenURL

  14. Garside H, Stevens A, Farrow S, Normand C, Houle B, Berry A, Maschera B, Ray D: Glucocorticoid ligands specify different interactions with NF-kappa B by allosteric effects on the glucocorticoid receptor DNA binding domain.

    Journal of Biological Chemistry 2004, 279(48):50050-50059. PubMed Abstract | Publisher Full Text OpenURL

  15. Ikeda A, Sun X, Li Y, Zhang Y, Eckner R, Doi T, Takahashi T, Obata Y, Yoshioka K, Yamamoto K: p300/CBP-dependent and -independent transcriptional interference between NF-kappa B RelA and p53.

    Biochemical and Biophysical Research Communications 2000, 272(2):375-379. PubMed Abstract | Publisher Full Text OpenURL

  16. Perkins ND: Integrating cell-signalling pathways with NF-κkB and IKK function.

    Nature Reviews Molecular Cell Biology 2007, 8:49-62. PubMed Abstract | Publisher Full Text OpenURL

  17. Salminen A, Ojala J, Huuskonen J, Kauppinen A, Suuronen T, Kaarniranta K: Interaction of aging-associated signaling cascades: Inhibition of NF-kappa B signaling by longevity factors FoxOs and SIRT1.

    Cellular And Molecular Life Sciences 2008, 65(7-8):1049-1058. PubMed Abstract | Publisher Full Text OpenURL

  18. Wadgaonkar R, Phelps K, Haque Z, Williams A, Silverman E, Collins T: CREB-binding protein is a nuclear integrator of nuclear factor-kappa B and p53 signaling.

    Journal of Biological Chemistry 1999, 274(4):1879-1882. PubMed Abstract | Publisher Full Text OpenURL

  19. Ashall L, Horton CA, Nelson DE, Paszek P, Harper CV, Sillitoe K, Ryan S, Spiller DG, Unitt JF, Broomhead DS, Kell DB, Rand DA, See V, White MRH: Pulsatile stimulation determines timing and specificity of NF-kappaB-dependent transcription.

    Science 2009, 324(5924):242-246. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. Fonslet J, Rud-Petersen K, Krishna S, Jensen MH: Pulses and Chaos: Dynamical Response in a Simple Genetic Oscillator.

    International Journal of Modern Physics B 2007, 21(23-24):4083-4090. Publisher Full Text OpenURL

  21. Krishna S, Jensen MH, Sneppen K: Minimal model of spiky oscillations in NF-κB signaling.

    PNAS 2006, 103(29):10840-10845. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  22. Nelson DE: Charaterisation and Manipulation of NF-κB Signal Oscillations in Single Living Cells. PhD thesis. School of Biological Sciences, University of Liverpool; 2005. OpenURL

  23. Matalka KZ, Tutunji MF, Abu-Baker M, Abu Baker Y: Measurement of protein cytokines in tissue extracts by enzyme-linked immunosorbent assays: Application to lipopolysaccharide-induced differential milieu of cytokines.

    Neuroendocrinology Letters 2005, 26(3):231-236. PubMed Abstract OpenURL

  24. Prabha C, Jalapathy KV, Matsa RP, Das SD: Role of TNF-alpha in host immune response in tuberculous pleuritis.

    Current Science 2003, 85(5):639-642. OpenURL

  25. Wang Y, Paszek P, Horton CA, Hong Y, White MRH, Kell DB, Muldoon MR, Broomhead DS: A systematic survey of the response of a model NF-κB signalling pathway to TNFα stimulation. [http://eprints.ma.man.ac.uk/1404/] webcite

    MIMS EPrint 2009.104 Manchester Institute for Mathematical Sciences, The University of Manchester; 2010. OpenURL

  26. Wiggins S: Introduction to Applied Nonlinear Dynamical Systems and Chaos. New York; London: Springer Verlag; 1990. OpenURL

  27. Pikovsky A, Rosenblum M, Kurths J: Synchronization : A Universal Concept in Nonlinear Sciences. Cambridge, UK: Cambridge University Press; 2003. OpenURL

  28. Nayfeh AH, Mook DT: Nonlinear Oscillators. New York: John Wiley & Sons; 1979. OpenURL

  29. Ermentrout B: Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students. Philadelphia: SIAM; 2002. OpenURL

  30. Doedel EJ, Paffenroth RC, Champneys AR, Fairgrieve TF, Kuznetsov YA, Oldeman BE, Sandstede B, Wang XJ: AUTO2000 : Software for continuation and bifurcation problems in ordinary differential equations. In Tech rep. California Institute of Technology, Department of Applied and Computational Mathematics, Pasadena CA 91125; 2000. OpenURL

  31. Arnold VI: Geometrical Methods in the Theory of Ordinary Differential Equations. 2nd edition. New York; London: Springer Verlag; 1988. OpenURL

  32. Arrowsmith DK, Place CM: An introduction to Dynamical Systems. Cambridge, UK: Cambridge University Press; 1990. OpenURL

  33. Denjoy A: Sur les courbes definies par les equations differentielles a la surface du tore.

    J Math 1932, 17(IV):333-375. OpenURL