### 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.

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:

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 *T*_{
R
} is a parameter indicating the presence or absence TNFα stimulation: *T*_{
R
} = 1 when the system is being stimulated with 10 ng/ml TNFα and *T*_{
R
} = 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 *T*_{
R
} = 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 *T*_{
R
} , allowing it to vary across the interval 0 ≤ *T*_{
R
} ≤ 1. In light of the modelling assumptions discussed above, one should not imagine that *T*_{
R
} depends linearly on the concentration of TNFα, but only that *T*_{
R
} 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 *T*_{
R
} as the parameter: Figure 3 is the resulting diagram. The model has a Hopf bifurcation (HB) at *T*_{
R
} = *T*_{*} ≈ 0.366, which means that for *T*_{
R
} > *T*_{
*
} the ratio of nuclear-to-cytoplasmic NF-*κ* B concentration exhibits sustained oscillations, but for *T*_{
R
} < *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.

### 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 (*T*_{
R
} = 1) stimulation alternate with periods during which *T*_{
R
} = 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.

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 < *T*_{
R
} < 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 \frac{1}{2} and \frac{3}{2} also appear. The first of these provides evidence that subharmonic resonance occurs over a wide range of pulsing frequencies.

When *T*_{
R
} ≥ 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 *T*_{
R
} = 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

{T}_{R}(t)=\epsilon (1+\eta \mathrm{sin}(2\pi \nu t)).

(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 *T*_{
R
} (*t*) has period τ = 1/*ν* and range

0\le \epsilon (1-\eta )\le {T}_{R}(t)\le \epsilon (1+\eta )\le 1.

The roles of the parameters are illustrated in Figure 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 *T*_{
R
} (*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 *T*_{
R
} (*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

where *p* and *q* are integers.

If *pν* and *qν*_{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

f=p\nu \approx q{\nu}_{0}

(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

which gives rise to the network of lines with slope of ±1. Finally, Figure 7 also exhibits sub-harmonic resonance: vertical strips for which *ν* ≈ *qν*_{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 *T*_{
R
} (*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 *ν* ≈ *mν*_{0}.

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}:

\tau =\frac{q}{\nu}\approx \frac{p}{{\nu}_{0}}

(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 *T*_{R} < 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 *T*_{
R
} < 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.