Modeling time delay in the NFκB signaling pathway following low dose IL-1 stimulation

Stimulation of human epithelial cells with IL-1 (10 ng/ml) + UVB radiation results in sustained NFκB activation caused by continuous IKKβ phosphorylation. We have recently published a strictly reduced ordinary differential equation model elucidating the involved mechanisms. Here, we compare model extensions for low IL-1 doses (0.5 ng/ml), where delayed IKKβ phosphorylation is observed. The extended model including a positive regulatory element, most likely auto-ubiquitination of TRAF6, reproduces the observed experimental data most convincingly. The extension is shown to be consistent with the original model and contains very sensitive processes which may serve as potential intervention targets.


Introduction
The transcription factor NFB is of central importance in inflammation and anti-apoptotic signaling. Upon stimulation of human epithelial cells with IL-1, NFB becomes activated due to proteasomal degradation of its cellular inhibitor IBα. This process requires phosphorylation of IBα by the upstream kinase IKKβ. Since sustained NFB-dependent expression of anti-apoptotic genes contributes to the maintenance of a range of cancers, its activity is tightly regulated and terminated by a negative feedback loop, as NFB promotes IBα synthesis. Accordingly, various approaches to anti-cancer strategies involve inhibition of the NFB signaling pathway [1].
Interestingly, NFB is converted into a pro-apoptotic factor upon stimulation with IL-1 + UVB. The persistence of this effect is ensured by sustained NFB activity [2] caused by sustained phosphorylation of IKKβ resulting in instant phosphorylation and proteasomal degradation of newly synthesized IBα. Chronic IKKβ phosphorylation, in turn, is due to UVB-induced inhibition of the responsible phosphatase PP2Ac [2]. We investigated the details of these processes using a systems biological approach, leading to the following ordinary differential equation model of IKKβ phosphorylation and dephosphorylation [3]: While IKKβ is very rapidly phosphorylated upon 10 ng/ml IL-1 stimulation, delayed IKKβ phosphorylation could be observed upon 0.5 ng/ml or 0.029 nM IL-1 stimulation (Figure 1). This slowly increasing phosphorylation activity is only insufficiently reproduced by the original model, so that a more detailed model of the signaling cascade appears appropriate for low IL-1 doses. However, the mechanism causing the signal delay is unknown to date. In principle, various reasons for signal delay are conceivable. Among the most prominent examples are double phosphorylation as occurring, e.g., in the MAPK cascade [4], mechanisms with irreversible inhibitors and positive feedback mechanisms. Here we will investigate these three mechanisms by extending the original model by each of them separately to predict the most probable explanation for the delay.

Modeling
Each of the three potential mechanisms delaying the signaling cascade following low doses of IL-1 is modeled as a separate building block (see Figure 2 The positive feedback mechanism describes the autocatalyzed activation of a protein T to the activated form Tu. This mechanism is given as In the double phosphorylation mechanism, the active form Ypp of a protein Y is reached after two sequential modifications (e.g., phosphorylations) catalyzed by ILRc: In the irreversible inhibitor mechanism, an irreversible inhibitor I blocks the activated form Xa of some protein X within the signaling cascade. The mechanism is implemented as While no specific irreversible inhibitor mechanism is known within the IL-1 signaling cascade, the double phosphorylation mechanism might, e.g., correspond to IRAK-1 or IRAK-4 phosphorylation [5]. A promising candidate for the positive feedback mechanism in the IL-1 signaling cascade is TRAF6: The E3-ligase TRAF6 is a key mediator of IL-1 induced signaling because its auto-ubiquitination is required to activate IKKβ [6]. Hence, the positive feedback mechanism might describe TRAF6 auto-ubiquitination. In the following, we therefore refer to this model extension as the TRAF extension. Note that, while no other comparable mechanism is known to date within the IL-1 signalling cascade, the TRAF module might also represent a different, yet unknown positively auto-regulated mechanism within the signalling cascade.
Experimental data for 0.5 ng/ml IL-1 stimulation were obtained from three western blot experiments measuring IKKβ Ser177/181 phosphorylation in human epithelial carcinoma cell line KB following IL-1 stimulation with and without UVB (300 J/m 2 ) costimulation until 120 min post-stimulation as described in [3]. They were subsequently extended with data from two western blots covering the first 30 min after IL-1 stimulation without UVB. Data for the 10 ng/ml IL-1 ± UVB stimulation are taken from [3]. Maximal IKKβ phosphorylation was clearly lower for low (0.5 ng/ml) than for high (10 ng/ ml) IL-1 stimulation, so that the maximal value could not be used as a reference value for scaling. Therefore, a separate scaling factor was used for IKKβ following 0.5 and 10 ng/ml IL-1 stimulation.
Modeling and analysis was performed with the Matlab (The MathWorks) based software tool PottersWheel (http://www.potterswheel.de freely available for academic use) [7], analogously to the procedures described in [3]. Particularly, the χ 2 value was used as objective function, with where M is the number of stimulations (i.e., M = 4), N i is the number of data points under stimulation i, y ij is data point j under stimulation i with standard deviation σ ij and y(t ij , θ) is the simulated value at time point j under stimulation i for the parameter vector θ. The upper bound of the scaling factor for low IL-1 stimulations was set to 4, as obtained from experimental data (not shown). Upper bounds for the remaining parameter values were derived as described in [3], lower bounds were not required.

Results and Discussion
Comparison of the different model extensions Fitting of the model with its different extensions to the experimental data yields better χ 2 values than the original model for all extensions ( Table 1). Comparison of the parameter values with the parameter values of the original model showed that the parameters downstream of the model extension, i.e., k uv , k p, and k dp , remain basically unchanged ( Table 2). While all model extensions are able to reproduce the delay (Figure 3; Figures S1 and S2 in Additional file 1), the TRAF extension is the model structure with the least degrees of freedom and lowest χ 2 value. Consequently, it also yields the best (lowest) value of the Corrected Akaike Information Criterion, AICc, which compares the goodness-of-fit of models with a different number of parameters (Table 1). In contrast, the other extensions only perform about as well as the original model in terms of AICc. The TRAF extension convincingly reproduces the experimental data and especially the signal delay ( Figure  3), and a biological counterpart to this mechanism exists within the IL-1 signaling cascade. In the following, we therefore focus on the TRAF extension. All parameter values of this extension are given in Table S1 in Additional file 1.

Detailed analysis of the TRAF6 model extension
Fixing of all model parameters also occurring in the original model to their value reported in [3] yields only slightly worse fit quality of χ 2 = 24.3. Remarkably, the fit quality is still better than the fit quality of the original model [3]   10 ng/ml IL-1 doses ±UVB are considered. As to the time course, the most prominent difference of the model extension compared to the original model consists in a longer plateau phase of the maximum concentrations, which is in better agreement to the measured data. In summary, the model extension is consistent to the model presented in [3] and justifies its additional complexity with a considerable enhancement of the fit quality and a consistent qualitative behavior for low doses. Although the parameter fixing performed above suggests that identifiability in the extended model might be worse than in the original one, identifiability analysis shows that all parameter values are identifiable from the given data: 2000 fits were performed with initial parameter values obtained by randomly perturbing the best parameter values by up to four orders of magnitude. The standard deviation of the parameter sets of the best 800 fits is not more than 1% for each parameter.
Furthermore, we calculated the relative sensitivity x against perturbations in parameter p i following IL-1 stimulation for two time-independent characteristics x Table 2 Values of the kinetic parameters downstream of the model extension in the original model [3] and in the extended models Model k p (s -1 ) k dp (s -1 ) k uv (s -1 ) Original ( of IKKp(t): the IKKp peak amplitude, x = max(IKKp(t)) and the mean of IKKp within the longest observed time period, For both sensitivity measures, a high dose dependency can be observed (Figure 4): for 10 ng/ml (0.588 nM) IL-1, parameter perturbations result in relatively weaker changes of the mean IKKp value (|s i | < 1), and have almost no effect on the peak amplitude. In contrast, the system is quite sensitive upon 0.5 ng/ml (0.029 nM) IL-1 stimulation: for several parameters, perturbations have an amplified effect on both peak amplitude and average IKKp value (|s i | > 1). Furthermore, the system is sensitive over a wide dose range, with sensitivities of up to ±5 (Figure 4).
Taken together, the positive feedback regulation presumably mediated via TRAF6 auto-ubiquitination [6] represents a consistent, well identifiable extension of the IKKβ phosphorylation model that reproduces the time delay at low IL-1 stimulations and is in accordance with literature findings.
On the other hand, it should be noted that determination of low levels of phosphorylated IKKβ is an experimentally challenging task, so that the values for 0.5 ng/ ml IL-1 stimulation should be regarded as somewhat semiquantitative measurements despite the partially low standard deviation. Against this background, though the postulated model consistently explains the experimental data one should neither completely reject the other two investigated mechanisms, which also perform well, nor lose sight of other potential delay mechanisms. Nevertheless, our results may provide an indication that the reported TRAF6 auto-ubiquitination [6] can indeed be interpreted as a positive feedback loop.
The very good accordance of the relatively small model with the experimental data together with the reported findings [6] make TRAF6 auto-ubiquitination a promising candidate for further research. Several system parameters could be characterized as very sensitive at low IL-1 doses. The respective processes therefore qualify as potential intervention targets in cancer therapy. The systems biological approach in combination with the necessary experimental validation can help to further elucidate important molecular mechanisms. Author details 1