A systematic survey of the response of a model NF-κB signalling pathway to TNFα stimulation.

White's lab established that strong, continuous stimulation with tumour necrosis factor-α (TNFα) can induce sustained oscillations in the subcellular localisation of the transcription factor nuclear factor κB (NF-κB). But the intensity of the TNFα signal varies substantially, from picomolar in the blood plasma of healthy organisms to nanomolar in diseased states. We report on a systematic survey using computational bifurcation theory to explore the relationship between the intensity of TNFα stimulation and the existence of sustained NF-κB oscillations. Using a deterministic model developed by Ashall et al. in 2009, we find that the system's responses to TNFα are characterised by a supercritical Hopf bifurcation point: above a critical intensity of TNFα the system exhibits sustained oscillations in NF-kB localisation. For TNFα below this critical value, damped oscillations are observed. This picture depends, however, on the values of the model's other parameters. When the values of certain reaction rates are altered the response of the signalling pathway to TNFα stimulation changes: in addition to the sustained oscillations induced by high-dose stimulation, a second oscillatory regime appears at much lower doses. Finally, we define scores to quantify the sensitivity of the dynamics of the system to variation in its parameters and use these scores to establish that the qualitative dynamics are most sensitive to the details of NF-κB mediated gene transcription.

Introduction

The transcription factor NF- is critical to the control of response to cellular stress and is involved in the regulation of cell-cycle/growth, apoptosis, inflammation and immunity (Hayden and Ghosh, 2008; Ghosh et al., 1998; Pahl, 1999; Gerondakis et al., 2006; Hoffmann and Baltimore, 2006; Hayden et al., 2006; Gerondakis et al., 2006; Bassères and Baldwin, 2006; Dutta et al., 2006). NF- is composed of homo- or hetero-dimers, with the ubiquitously expressed RelA:p50 hetero-dimer being the primary inflammatory mediator (Hoffmann and Baltimore, 2006).

In the absence of any stimulus, NF- is held in an inactive state in the cytoplasm where is sequestered by association with inhibitory () proteins including , , and . In response to stimulation by cytokines, including , activated kinases (IKKs) phosphorylate the proteins, targeting them for degradation via the ubiquitin-proteasome pathway (Ghosh et al., 1998). Liberated NF- translocates to the nucleus and regulates target gene transcription, including highly inducible and zinc finger protein A20 genes (Scott et al., 1993; Song et al., 1996). This transcriptional control constitutes negative feedback regulation (Lipniacki et al., 2004). Newly synthesised binds to nuclear NF-, leading to export of the complex to the cytoplasm (Arenzana-Seisdedos et al., 1995), while A20 inhibits the NF- signalling cascade (Lipniacki et al., 2004; Wertz et al., 2004), acting upstream of IKK.

Coupled negative feedback loops may lead to oscillations, and, indeed, Hoffmann et al. (2002) used electromobility shift assay (EMSA) to observe damped oscillations in NF- nuclear activity at the population level. Subsequently White's lab (Nelson et al., 2004) observed sustained oscillations in NF- nuclear localisation at the single-cell level using fluorescence microscopy. In these experiments cells were stimulated continuously with a high dose of (10 ng/ml). In contrast, blood plasma measurements suggest that physiological concentrations of may be considerably lower (Nakai et al., 1999; Damas et al., 1989; Matalka et al., 2005). Recent single-cell data from White's and other labs demonstrated that a fraction of cells in the population can respond to concentrations of as low as few pg/ml, albeit with an apparently stochastic delay (Turner et al., 2010; Tay et al., 2010; Kingeter et al., 2010).

Dynamical responses of the NF- system regulate many physiological processes associated with inflammatory signalling. Frequency of oscillations, their persistence and other characteristics are fundamental in controlling patterns of downstream gene expression (Ashall et al., 2009; Tay et al., 2010). Therefore, a quantitative understanding of systems responses at physiological doses of the stimulus is required to elucidate the biological function of these dynamics. Here we use numerical bifurcation theory to survey systematically the dynamical responses of a model of the NF- system developed by Ashall et al. (2009). First, in Section 2 we consider the one parameter bifurcation problem associated with variation in the intensity of stimulation. Of course, the dynamics of the system would also be expected to change if one perturbed, for example, the reaction rates. In order to understand the influence of the variation of these rates on the qualitative dynamics of the system, we also carried out an extensive programme of two-parameter bifurcation studies that are discussed in Section 3. Finally, in Section 4, we define a parameter sensitivity score to quantify the sensitivity of the dynamics to variations in each parameter and use these scores to inform a discussion of further model development.

One-parameter bifurcations

The deterministic model in Ashall et al. (2009) is among the more detailed models of the NF- signalling system and includes, for example, explicit terms for mRNA concentration. It stands in contrast to reduced models such as those of Krishna et al. (2006) and Fonslet et al. (2007) that aim to capture essential mechanisms at the expense of some biochemical detail. It consists of two negative feedback loops: one in which NF- transcriptionally regulates its inhibitor and another in which it regulates expression of A20 protein, which acts indirectly to suppress the kinase IKK involved in liberating NF- from –NF- complexes. The model network is illustrated in Fig. 1 and specified in detail by a system of ordinary differential equations (ODEs) listed in Appendix A along with the associated parameters, which are listed in Table 1.

Fig. 1

A graphical representation of the model NF- signalling network detailed in Appendix A. Here ovals represent proteins and the two rounded rectangles represent messenger RNAs whose production is regulated by NF-.

Table 1

Ranges of model parameters and percentages-of-range for which the one-parameter bifurcation diagram in T has the same topology as that in Fig. 2.

Parameter		Range surveyed	Value in model	% one HB-point
Num.	Name
1	TR	(0, 1)	1.0	—
2	totIKK	(0.001943, 0.21102)	0.08	97.6
3	totNFkB	(0.032147, 0.1 )	0.08	59.8
4	kv	(1.7658, 5.7363)	3.3	46.5
5	kp	(1.00E−05, 0.0015968)	0.0006	98.3
6	ka	(0.00010279, 2.9586)	0.004	97.4
7	ki	(0.0007336, 0.11817)	0.003	75.5
8	kda	(9.48E−05, 0.021019)	0.0045	97.9
9	k_bA20	(6.81E−06, 0.0093601)	0.0018	99.6
10	kc2	(0.00028691, 0.46452)	0.074	99.6
11	kc3	(0.0029809, 1.07)	0.37	99.2
12	kdegpin	(0, 0.736)	0.1	100.0
13	ka1	(0.10325, 1.629)	0.5	79.4
14	kd1	(0.0002, 0.01128)	0.0005	60.0
15	kdegf	(0.0004, 0.0014042)	0.0005	10.0
16	kdegc	(0, 0.00019312)	0.000022	100.0
17	ki1	(0.00155, 0.047406)	0.0026	40.4
18	k_e1f	(0, 0.00021787)	52E−6	100.0
19	k_e1c	(4.48E−05, 0.38702)	0.01	99.6
20	ki2	(6.84E−05, 0.0015916)	0.00067	89.8
21	ke2	(0, 0.016038)	3.35E−4	100.0
22	h	(0.68071, 2.4335)	2	21.7
23	k	(0.055, 0.11369)	0.065	15.4
24	kitria	(7.06E−08, 1.68E−07)	1.4E−7	19.8
25	ktria	(0.20438, 0.6)	0.5	20.0
26	kdegtia	(4.98E−06, 0.00059088)	0.0003	97.0
27	kitra	(5.13E−08, 2.63E−06)	1.4E−7	63.4
28	ktra	(0.18323, 11.436)	0.5	63.4
29	kdegta	(3.30E−05, 0.0020154)	0.00048	93.1

This model assumes that IKK exists in one of three forms: neutral denoted by IKK (the form that can be activated), active denoted by IKK (the form that phosphorylates and its complexes) and inactive denoted by IKK (a state that cannot be activated). The model also assumes that inactive IKK can be recycled into active IKK and that A20 acts by inhibiting this recycling process.

The full model contains 14 ordinary differential equations, including equations governing the concentrations of the three forms of IKK; the free forms of , NF- and A20; the –NF- complex and the phosphorylated forms of and –NF-; and A20 and transcript. Most biochemical reactions are modelled with mass-action kinetics: the exceptions are the NF- transcriptional regulation of the A20 and genes and the action of A20 on IKK recycling. The regulation of gene expression by NF- is modelled by an increasing Hill function of the formwhere x represents the concentration of the transcription factor, k is the concentration required to produce a half-maximal response, and gives the maximal rate of transcription. The parameter h is often called the Hill coefficient and it regulates the nonlinearity of the response. The graph of f(x) in (1) is an increasing S-shaped curve that becomes more step-like as one increases h. The inhibition of IKK to IKK cycling mediated by A20 is modelled by a decreasing Hill function

Although the model describes the temporal evolution of 14 chemical concentrations, there are really only 12 independent quantities as the system has two conservation laws. First, NF- is neither created nor destroyed, so the total amount of NF- (comprising free NF- and its complexes, both cytoplasmic and nuclear) remains constant. Similarly, total IKK (encompassing the neutral, active and inactive forms) is also conserved. Additionally, the concentration of phosphorylated (pIkB in the notation of Ashall et al., 2009) does not appear in any of the ODEs except the one involving its own degradation. Thus we can ignore it and, taking account of the two conservation laws, reduce the model to a system of 11 ODEs: these are given in Appendix A.

The bifurcation and sensitivity analyses reported below are based on this reduced system, in which the intensity of stimulation, denoted by T, and the total amounts of NF- and IKK are treated as parameters. The latter two are denoted by totalNFkB and totalIKK, respectively. The reduced model – when used with the parameters in Table 1 – captures the two main features observed in the single cell experimental data in Ashall et al. (2009) and Nelson et al. (2004): (i) continuous, high dose (T=1) stimulation leads to sustained oscillations with a period near 100 min and (ii) in the absence of stimulation, the system remains in a steady state and does not oscillate.

Varying the dose

In the model, activates NF- via activation of IKK. Eqs. (3) describe the rates of the neutral and active IKKwhere are parameters, is the concentration of the protein A20 and, in Ashall et al. (2009), T indicates the presence or absence of high dose (10 ng/ml) stimulation: T=1 when is present and T=0 when it is absent. It is thus natural to generalise the role of T and to use values to represent the range of intensities of stimulation. Note that we do not assume that the value of T depends linearly on the concentration of , but only that the relationship is monotone and saturates at T=1 for high dose stimulation.

Now consider T as a bifurcation parameter. The parameter totalNFkB corresponds to the initial concentration of complex and is set to , while the concentrations of all other forms of NF- are set to zero. Similarly, totalIKK is set by assigning the initial concentration of neutral IKK to and the concentrations of all other forms of IKK to zero. The remaining parameters are held at the values used in Ashall et al. (2009).

Starting from the steady state equilibrium associated with unstimulated cells (T=0), the numerical bifurcation using XPPAUT indicates a branch of steady state solutions that changes stability by way of a supercritical Hopf bifurcation (HB) when . This means that, with other parameters held fixed, if any initial configuration of concentrations will undergo damped oscillations and eventually settle into a steady state with constant concentrations. But if , the system will evolve into sustained oscillations. This is illustrated in the left panel of Fig. 2, while the right panel shows the T dependence of the period of oscillation.

Fig. 2

Left panel: bifurcation diagram for T, where the thick black curve represents a branch of steady states, the thin black curve represents a branch of unstable steady states and the red curve represents a branch of limit cycles. The upper part of the red curve shows the peak values of the oscillations while the lower part shows the troughs. Right panel: period of the limit cycle as a function of T.

Sustained oscillations eventually disappear via a second Hopf bifurcation at T=156, but such a large value of T is not biologically meaningful. We emphasise again that T is not proportional to concentration: rather, it indicates the extent to which the signalling machinery is saturated. Single-cell measurements show that the system is already saturated at doses of 10 ng/ml (Tay et al., 2010; Turner et al., 2010) so we have restricted our analysis to values of T between 0 and 1.

From the right panel of Fig. 2 we see that the period of the limit cycle varies over a small range, increasing from 93 min to 99 min as T increases from 0.366 to 1. We note that when , experimental observations might still appear to indicate the presence of oscillations, as all the stable steady states in the branch are foci (the Jacobian along the branch of steady state has complex eigenvalues). For levels of stimulation in this particular region, the oscillations will damp away exponentially, but – especially for T near the HB-point – the damping rate is very small, so that oscillations in this regime, though decaying slowly, may appear to persist throughout the period of observation. Additionally, molecular noise can support sustained stochastic oscillations even for T below the HB point: see Turner et al. (2010).

Two-parameter bifurcations

As the values of parameters in the model are determined either by direct measurement or via model-fitting, both of which are subject to error, it is natural to ask how the structure of the one-parameter bifurcation is influenced by the changes to the other parameters. One can begin to address this question by performing two-parameter bifurcation analyses. Consider, for example, the parameter totalIKK: in Section 2.1 we obtained a one-parameter bifurcation diagram for the case . Now if we change value of totalIKK by a small amount and recompute the one-parameter bifurcation diagram, then typically we will get a new diagram similar to Fig. 2, but with the HB-point in a slightly different position. This suggests that the one-parameter bifurcation diagram could be summarised in terms of the existence and location of the HB-point. Further, the dependence of the HB-point's position on totalIKK can be summarised by a two-parameter bifurcation diagram, which is a projection of all bifurcation points into the (T, totalIKK)-plane. Such a diagram shows how the bifurcation point moves with variation in a second parameter value, which here is the value of totalIKK.

Two-parameter bifurcation diagrams can also be computed with the numerical bifurcation package AUTO (Doedel et al., 2000): Fig. 3(left) is the two-parameter diagram for T and totalIKK. In the left panel we can see that near the original HB-point, a cell with a higher level of totalIKK would require stronger stimulation to support sustained oscillations. But the slope reverses for : there a cell with a lower concentration of total IKK needs stronger intensity of stimulation to support oscillations. And if , then no sustained oscillations occur for any level of stimulation.

Fig. 3

The two-parameter bifurcation diagrams for T and totalIKK (left panel) and T and totalNFkB (right panel). Here each point on the curves represents a HB-point and the curves divide the parameter spaces into two types of region: there is a limit cycle for each pair of the parameters in Region 1, in which and there is a non-oscillatory steady-state for each pair of parameters in Region 2.

One might expect that if the perturbation to the second parameter is large, novel dynamical behaviour may appear. This happens when we increase the value of total NF- from to a value greater than : see the right panel in Fig. 3. When, for example, , the corresponding one-parameter bifurcation diagram, which appears in the left panel of Fig. 4, has three HB-points.

Fig. 4

The bifurcation diagram for T when : there are three HB-points at , 0.01 and 0.366, respectively. Sustained oscillations occur for values of T in either of the intervals (0.002, 0.01) and (0.366, 1]. The right panel shows the period of these oscillations.

The system has sustained oscillations when T lies in either of the intervals or . The amplitude of the oscillations is quite small when T is in the interval with a period in the range of .

Bifurcation structures similar to those in Fig. 4 – with three Hopf bifurcations – also appear in the two-parameter diagrams for T and each of the parameters h, totalNFkB, k, k, k, k, , k, k and k. All remaining two-parameter diagrams appear in Appendix C, while Table 1 lists the ranges of parameter values over which the one-parameter bifurcation diagram in T retains the same structure as is illustrated in Fig. 2—that is, a single supercritical Hopf bifurcation and no others.

The bifurcation diagrams in Appendix C make it clear that a wide range of parametric variation produces only two types of two-parameter bifurcation diagrams. One type topologically resembles the diagram on the left panel of Fig. 3 where, for a given value of the second parameter, there is a unique critical value of T separating the levels for stimulation that induce sustained oscillations for those that produce only damped oscillation. The second type of bifurcation diagram resembles that in Fig. 3: variation in of second parameter can lead to oscillating behaviour associated with the high levels of stimulation, but for suitable parameter values a second, separate oscillatory regime exists for very low-intensity stimulation. These low-dose oscillations are isolated from the high-dose oscillations in the sense that there are intermediate values of T for which no oscillations occur. Further, the low-dose oscillations are readily distinguishable from those in the adjacent damped regime: see Fig. B1 in Appendix B.

Fig. B1

A typical experiment tracks single cells for  min, which is long enough to distinguish the low-dose oscillatory regime (blue curve) from the exponentially damped oscillations that occur at slightly higher (black curve) and lower (green curve) levels of stimulation. The model here is the same as that depicted in Fig. 4 and has . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Parameter sensitivity analysis

In this section we define a sensitivity score to quantify the sensitivity of the model to changes in the parameters. In the previous section we observed that for parameters in the ranges listed in Table 1, the one-parameter bifurcation diagram in T retains the qualitative structure illustrated in Fig. 2: there is a single HB-point. Our sensitivity score measures the impact of variation in the other parameters on the location of the HB-point in the associated one-dimensional T diagram. In addition, as changes in the period of the oscillation may lead to changes in the pattern gene expression (Ashall et al., 2009), we also investigate the sensitivity of the period to variation of the parameters.

Consider a parameter (other than T) whose unperturbed value is . The associated sensitivity score is defined as the average value of with in a region where is a prescribed fractional change in the parameter and is the value of T at the Hopf bifurcation point as described in Section 2.1.

We take the following steps to compute the score:

Read the pairs of from the two-parameter bifurcation diagram over the region : this bifurcation curve and the corresponding pairs comes from AUTO.

Use piecewise cubic Hermite interpolation to fit a continuous function to the pairs.

Sample uniformly-spaced points from the interpolated curve, where the number of sampling points m is the same for all parameters. Denote these points as: , with for .

The sensitivity score is then defined by the approximate integralwhere is the parameter value used in Ashall et al. (2009) and is the critical value from Fig. 2.

We computed sensitivity scores by varying each parameter individually by a fraction of : they are displayed in Fig. 5, where the numbering of the parameters is the same as that used in Table 1. We found that the parts of the two-parameter bifurcation curves spanned by these 10% variations were well-approximated by line segments, so based our scores on samples of m=21 points for each parameter, which is more than sufficient to guarantee convergence of the integrals.

Fig. 5

Sensitivity scores for the 26 reaction rates and the two concentration parameters (totalNFkB and totalIKK) constructed using Eq. (4) with : (left) absolute value of the sensitivity score against parameter number (see Table 1) and (right) sensitivity score against the index of parameters.

The scores are clearly divided into two groups: those with higher scores and those with lower ones. Using a threshold value of , parameters with the higher sensitivity scores are listed in Table 2, which shows that the system is most sensitive to variation of the parameters related to transcription and translation of and A20 and to those influencing IKK activity. These results are consistent with those obtained earlier by Yue et al. (2006, 2007, 2008) and Ihekwaba et al. (2004). The analysis in these earlier papers was based on a deterministic model developed by Hoffmann et al. (2002), which includes only the core NF- and feedback loops. While Yue and Ihekwaba analysed the sensitivity of a certain transient response to the variations in the model's parameters, we have here used both a more recent model and have studied the sensitivity of long-time dynamical behaviour—the question of whether there are stable oscillations or not. Our results thus support and complement the earlier work, bolstering their finding about the high sensitivity of the dynamics to variation in those parameters related to transcription and translation of as well as those influencing IKK activity. Further, our analysis indicates a role for the parameters associated with the transcription and translation of A20 that is comparable in importance to those influencing production and dynamics of .

Table 2

The parameters with higher sensitivity scores, arranged in order of decreasing sensitivity.

Parameter		Description	Score
Num.	Name
24	k	Nuclear NF-κB concentration at half-maximal transcription	1.17
4	kv	Ratio of cytoplasmic to nuclear volume	−1.01
23	h	Order of Hill function controlling transcription	−0.987
25	kitria	Inducible IκBα mRNA synthesis	−0.984
26	(k_tria)	IκBα Translation rate	−0.984
2	totIKK	total IKK	0.914
7	ki	Spontaneous IKK inactivation	−0.912
27	kdegtia	IkBa mRNA degradation	0.91
5	kp	IKKn production rate	0.87
8	kda	Constitutive A20 degradation	0.87
28	kitra	Inducible A20 mRNA synthesis	−0.86
29	ktra	A20 mRNA translation rate	−0.86
9	k_bA20	A20 inactivation of IKK  (All TNF conditions)	0.85
11	k_c_3	Catalysis of IKK.IkBa.NFkB trimer	0.79

Since the sensitivity score defined here has the same sign as the mean values of the slopes along the two-parameter bifurcation curves, and as the relevant segments of the curves are nearly linear, the sign of the scores provides further information. For example, note that k, k, k and k all have negative scores. This implies that the faster the transcription or translation of and A20, the lower the dose of needed to induce sustained oscillations. Similar properties can be read off for other parameters from Fig. 5.

As pointed out in Ashall et al. (2009), the period of the oscillations plays an important role in selecting which of the genes that NF- targets are expressed. We thus also considered the sensitivity of the period of the NF- oscillations (while T=1) to the changes in parameters. That is, fixing T=1 and varying each of the other parameters by 10%, we compute the sensitivity scores by substituting T by the value of the period at the corresponding parameters in (4). Using the same sampling principles, we obtain the resulting scores that are shown in Fig. 6.

Fig. 6

Scores measuring the sensitivity of the period of the NF- oscillations to the changes in parameters: (left) absolute value of sensitivity and (right) sensitivity score against parameter index.

Note that all scores in Fig. 6 are much smaller than 0.5, which is the threshold we used for distinguishing the higher and lower score groups in Fig. 5. This that, compared to the location of the Hopf bifurcation, the oscillatory period is much less sensitive to the parameters.

Discussion

The transcription factor NF- is pivotal in controlling body's innate immunity, but can also contribute to carcinogenesis (Karin et al., 2002). Live-cell imaging has established that NF- undergoes sustained oscillations in cellular localisation that can be induced by continuous, high doses of stimulation (Ashall et al., 2009; Nelson et al., 2004). More recently, the remarkable ability of the NF- system to oscillate in response to much lower, physiological doses of has also been demonstrated experimentally (Turner et al., 2010; Tay et al., 2010; Kingeter et al., 2010). The data from White's lab (Turner et al., 2010) showed that individual SK-N-AS neuroblastoma cells may oscillate in response to stimulation with doses as low as picomolar , though this response is apparently probabilistic, with the fraction of responsive cells declining with dose.

Other studies using different cell types showed similar stochastic activation of the pathway. However, oscillations of the transcription factor at lower doses were less apparent (Tay et al., 2010; Kingeter et al., 2010). Stochastic models of the system suggest that probabilistic nature of low-dose responses may arise from the noisy activation of the transduction pathway leading from the receptor to IKK module (Turner et al., 2010) and more specifically, to the limited availability of trimers for receptor binding (Lipniacki et al., 2007; Tay et al., 2010). Because of their stochastic nature, none of these models is straightforwardly amenable to the sort of systematic bifurcation analyses reported above.

Here, we have surveyed the qualitative behaviour of a recent deterministic model of the NF- system (Ashall et al., 2009), and used bifurcation analysis to characterise all possible responses to stimulation with the intensity in the range of 0–10 ng/ml. Further, we have quantified the impact of variation in the remaining parameters on the level of stimulation required to induce oscillations. Using the parameter T – a proxy for the intensity of stimulation – as bifurcation parameter, we found that sustained oscillations appear via a supercritical Hopf bifurcation. That is, there exists a critical intensity of stimulation such that sustained oscillations will occur whenever , while only damped oscillations will occur for . Further, a similar bifurcation structure persists across a wide range of variation in most of the model's other parameters. There is, however, a group of parameters – about a third of those appearing in the model – for which sufficiently large modulations introduce the possibility of a second, distinct sort of periodic oscillation appearing at very low levels of stimulation.

Building on our bifurcation analyses, we performed sensitivity analyses to quantify the effect of parameter variation on the location of the supercritical Hopf bifurcation and on the period of the oscillation. The location of the bifurcation proved highly sensitive to variation in those parameters associated with transcription and translation of and A20 as well as those associated with IKK activity, a result consistent with previous work (Ihekwaba et al., 2004; Yue et al., 2006, 2007, 2008). By contrast, the period was much less sensitive to parametric variation.

The low-dose limit-cycle regime described above appears for , which corresponds to levels of stimulation that induce less than 1% of a saturating response. An oscillatory response to such feeble stimulation offers a tantalising suggestion of a connection to data from low-dose experiments, but there are a number of obstacles to making this connection concrete. As we emphasised in Section 2, our main bifurcation parameter T is only a proxy for the intensity of stimulation: it is an abstraction representing the consequences of that part of the signalling cascade that extends from the receptors at the cell's surface down to IKK. A modest extension to the model might treat the upstream part of the signal transduction chain as a Hill function (1) whose input would be extracellular and whose output would play the role of T. This suggests an explanation for the tenfold range of T variation over which the current model predicts no oscillations at all: these values could correspond to the step in the Hill function and thus to a narrow range of concentrations that might prove difficult to locate experimentally. But as neither we nor our collaborators have data that would allow us to fit such an extended model, we have not pursued it any further.

The emergence of the low-dose oscillatory regime also depends on modulation of some of the model's other parameters, including expression levels of NF- and IKK, rates of transcription and transport among others (Figs. 3, C1 and C2). Data available in the literature show that during carcinogenesis many components of the NF- signalling systems are changed or mutated in such a way that the activity of the transcription factor is increased (for example in multiple myeloma, Annunziata et al., 2007; lung cancer, Tang et al., 2006; and colon cancer, Charalambous et al., 2009). These measurements are not really quantitative, and thus direct comparison with our model is not possible, but one could perhaps speculate whether during cancer progression parameters of the NF- system change in such a way as to promote anomalous limit cycle oscillations in a low-dose tissue micro-environment. This could enable more efficient proliferation of cells, perhaps by driving oscillation-dependent interactions with a cell cycle system – for example, via cyclin D (Sée et al., 2004) – that are inhibited or not activated in healthy tissue.

Methods

Bifurcation diagrams were prepared with XPPAUT (Ermentrout, 2002) and the command line interface Rauto (Schilder, 2007): both rest on the powerful numerical bifurcation package AUTO (Doedel et al., 2000). Sensitivity scores were computed with MATLAB R2006b: code is available on request.