Computational Models of the NF-κB Signalling Pathway

In this review article, we discuss the current state of computational modelling of the nuclear factor-kappa B (NF-κB) signalling pathway. NF-κB is a transcription factor, which is ubiquitous within cells and controls a number of immune responses, including inflammation and apoptosis. The NF-κB signalling pathway is tightly regulated, commencing with activation at the cell membrane, signal transduction through various components within the cytoplasm, translocation of NF-κB into the nucleus and, finally, the transcription of various genes relating to the innate and adaptive immune responses. There have been a number of computational (mathematical) models developed of the signalling pathway over the past decade. This review describes how these approaches have helped advance our understanding of NF-κB control.


Introduction
Cells use complex networks of interacting molecular components to transfer and process information. It is increasingly appreciated that the robustness of various cellular processes is rooted in the dynamic interactions among the cell's many constituents [1], such as proteins, DNA, RNA and small molecules. Barkai and Leibler [2] advise that the "computational devices of living cells" are responsible for many important cellular processes, including metabolic pathways, signalling pathways and regulatory mechanisms.
Intracellular pathways possess a high degree of complexity, with the role of their constituent components not only being defined by their individual function, but also by their interaction with other components [3]. It has recently emerged that many biological components are shared by several pathways, which merge into large and complex interconnected networks when viewed from a systems perspective. Oltvai and Barabasi [4] propose that we build on this approach and believe that by viewing a cell as an interconnected network of genes and proteins, we will be able to develop viable strategies for addressing the complexity of living systems. The field of systems biology is concerned with understanding this complexity by focusing on the interactions between the various components of living cells [5]. It has more recently become associated with the use of an integrative approach, combining experimental and computational modelling to identify the molecular mechanisms underlying the properties and behaviours of complex biological systems, termed computational biology [6,7].
The discipline of computational biology encompasses all fields of biology that have joined forces with computer science and engineering. Just as biology can be broken down into sub-disciplines, such as cell biology, biochemistry, immunology, etc., a growing specialty within the discipline of computational biology relates to the immune system and has been termed computational immunology. The immune system provides our primary protection mechanism against invading organisms (e.g., bacteria, viruses, fungi or other microorganisms) and cells that have become uncontrollably changed [8]. Many of the responses following the identification of a foreign organism or toxic molecule are destructive in nature and facilitate the containment or removal of the pathogen [9]. Following recognition of a pathogen, a plethora of immune responses develop over time, which are caused by modification to the activity of intracellular enzymes or activation factors, e.g., transcription factors, regulating expression of genes involved in host defence. There are a number of well-defined signalling pathways in humans that lead to activation of gene expression. Of these, the nuclear factor-kappa B (NF-κB) pathway is believed to be central to the regulation of inflammatory responses, as scientists have not yet discovered an inflammatory gene that is not controlled, at least in part, by NF-κB.
The aim of the present review is to discuss the development of computational modelling of the NF-κB signalling pathway and how this has facilitated significant advancements in our understanding of its regulation.
NF-κB results in the induction of gene expression of a number of inflammatory response proteins, including the inhibitor IκBα, thus providing a negative feedback mechanism that controls transient activation of the pathway. Figure 1. High-level representation of how the multitude of inducers of the NF-κB signalling pathway, which converge on the NF-κB signalling module (NF-κB, IκBα and IKK), before diverging to regulate a multitude of gene products, has the topology of a bow-tie motif from systems biology. Pahl [11] advises that there are over 150 known inducers of NF-κB activation, which can be allocated to a number of categories of environmental stimuli. Similarly, Pahl advises that more than 150 genes are known to be regulated by NF-κB, which can again be categorised into a number of key physiological responses.
In addition to the various NF-κB dimers, IκB inhibitor proteins and IKK complexes, there are also a large number of upstream components that facilitate signal transduction through the pathway. Of these, the proteins within the cell receptor complex are the pivotal components in facilitating the cells response to its environment. Most surface receptors are transmembrane in nature and contain an extracellular, ligand-binding domain, responsible for sensing the incoming signal, and an intracellular domain, responsible for initiating the propagation of the signal through the cascade of different kinases until it reaches the IKK-NF-κB-IκB signalling module. Along with the cell membrane receptor proteins, the cell receptor complex may also contain co-receptors and adaptor proteins. For example, following ligation with lipopolysaccharide (LPS), the Toll-like receptor 4 (TLR4) protein dimerises and associates with CD14 (a co-receptor for LPS) and an additional cellular protein MD2 [31] in order to signal the presence of the bacterial product. Furthermore, some proteins within a receptor complex mediate the signal transduction between proteins by bringing them together and are referred to as adaptor proteins. They function as docking proteins to co-locate signalling molecules into multiprotein signalling complexes. The formation of the receptor-co-receptor-adaptor protein complex and related events at the level of the receptor have a primary role in signal transduction. Along with the core components of the signalling network discussed above, another set of components involved in regulating gene expression are a group of proteins termed modulators. Modulators are used within signal transduction networks to modify the activity of transcription factors and often act by binding to them directly and, thus, affect the expression of target genes [32]. With particular reference to the NF-κB signalling pathway, Li et al. [33] have recently used computational modelling and analysis of gene expression profiles to predict a number of key modulators and the associated target genes of NF-κB. They identified 365 candidate modulators of the RelA monomer, including 334 that had not been previously reported. They developed a probabilistic model of the regulatory network and discovered a number of functional subnetworks, which were believed to contain new pathways for modulating RelA transcriptional activity and specificity. Figure 3. Unified Modelling Language (UML) class association diagram for the canonical NF-κB signalling pathway. The layout used follows the spatial aspects of the cascade of component relationships involved in the biological system from cell membrane activation down to gene transcription within the nucleus. The diagram conveys directed relationships by the labels at each end of the connectors. The arrowheads pointing towards the nuclear receptor convey inheritance, such that exporting receptor and importing receptor are sub-types. The unfilled diamonds on the NF-κB-IκB complex denote aggregation, i.e., an IκB component and an NF-κB component aggregate to form the inhibited complex.

Existing Computational Models of NF-κB
A significant body of knowledge has been generated through wet-lab experimentation since the discovery of NF-κB in 1986. Following the move by a large section of the scientific community away from the reductionist approaches and towards a systems-level understanding, the use of computational modelling and simulation has been adopted by a growing number of research groups investigating NF-κB. Because of the central role of the NF-κB inhibitor in the control of transcriptional activation, its regulation has been the focus of the majority of in silico models of pathway regulation. Recently, equation-based and agent-based models have been used within a predictive capacity to generate hypotheses for testing through additional wet-lab experimentation.
We will review existing computational models that have provided the major advancements throughout this section. Figure 3 depicts a UML (Unified Modelling Language) class association diagram for the canonical signalling pathway, which has formed the basis for the majority of the computational models to date.

Deterministic Differential Equation Models
The first mathematical model of NF-κB dynamics was developed by Carlotti et al. [34]. This model used reaction kinetics, built with Microsoft Excel and Visual Basic, to focus on IκBα association and dissociation rates, along with IκBα and NF-κB translocation between the cytoplasm and the nucleus. The first model of the wider NF-κB signalling pathway, however, was that of Hoffmann et al. [35]. Hoffmann et al. developed an ODE-based computational model using Gepasi software (version 3.1; Mendes [36,37]). They were particularly interested in the temporal control of NF-κB activation by the coordinated degradation and synthesis of IκB proteins. The model therefore focuses on reactions that govern IκB dynamics, including synthesis and degradation, along with cellular localisation and association/dissociation with NF-κB.
The Hoffmann model comprised a system of 24 ODEs describing reaction kinetics equations of the change in concentration (with respect to time) of cytoplasmic and nuclear NF-κB, cytoplasmic and nuclear IκB (-α, -β and -), cytoplasmic IKK and resultant complexes that they may form. The model also contained 30 parameters, which were gained using a combination of biochemical experimentation within their lab, a review of published literature and estimation through model fitting techniques. The model describes the differential functions of IκBα, β and isoforms on the regulation of NF-κB activation and used knockout mice for the three isoforms, cross-bred to yield double-knockout cells where NF-κB was inhibited by a single IκB isoform. These cells were stimulated with TNFα in order to induce the NF-κB signalling response. Their conclusions from wet-lab experimentation were that coordinated degradation, synthesis and localisation of all three IκB isoforms is required to generate the characteristic NF-κB activation profile.
Three computational models were developed to represent these double knockouts: β -/--/-(IκBα), α -/--/-(IκBβ) and α -/β -/-(IκB ), which were calibrated so that simulation results were consistent with wet-lab dynamics. Simulation yielded two very different results, predicting that: (i) IκBα mediates rapid NF-κB activation and strong negative feedback regulation, resulting in an oscillatory NF-κB activation profile; and (ii) IκBβ and IκB respond more slowly to IKK activation and act to dampen the long-term oscillations of the NF-κB response.
The authors subsequently built on this initial work by performing new wet-lab knockout experiments to isolate cleanly the endogenous free-and bound-IκB protein pools and probed their degradation with kinase knockouts and pharmacological inhibitors [38]. They used the previous model to investigate the steady-state regulation of the NF-κB signalling module and its impact on stimulus responsiveness, through manipulating rate constants of IκB isoform degradation through either IKK-mediated or IKK-independent processes. This in silico experimentation revealed a homeostatic NF-κB signalling module in which differential degradation rates of free and bound pools of IκB represent a novel cross-regulation mechanism that imparts functional robustness to the signalling module.
Recently, Cheong et al. [39] have reviewed this initial work and argue that the functions of the three IκB isoforms combine to allow the signalling module to distinguish between short and long lasting stimuli. Furthermore, Kearns et al. [40] through reimplementation in MATLAB and augmentation of the model, discovered that IκB provides negative feedback to control NF-κB oscillatory dynamics, but that a delay in IκB transcription renders this negative feedback in antiphase to that of IκBα. The authors therefore suggest that IκB has a role in dampening the IκBα-mediated oscillations during long-lasting NF-κB activity. This must be tempered, however, as subsequent work by Paszek et al. [41] (see the semi-stochastic modelling section) suggested that the antiphase negative feedback by IκBα and IκB is optimised to generate heterogeneous NF-κB oscillations between cells and, thus, a population-level robustness to external perturbations.
Following additional wet-lab experimentation, they also hypothesise that the relative strength of IκBα and IκB feedback mechanisms and their temporal relationships to each other may account for cell type-specific regulation of NF-κB dynamics. Basak et al. [42] have also updated the model to introduce the p100 protein, which is part of the non-canonical signalling pathway, and LPS induction of IKK-mediated IκB degradation. They showed that p100 acts as an inhibitor of NF-κB activation and termed this IκBδ, as it was found to be a bona fide IκB protein. They used the model to explore feedback regulation and dynamics of p100 and suggest that there may be crosstalk between the canonical and non-canonical signalling pathways. More recently, Shih et al. [43] have further augmented the model of Basak et al. and found that the newly discovered inhibitor IκBδ provides negative feedback in the presence of persistent, pathogen-triggered signals, by dampening NF-κB responses during sequential stimulation events. The authors subsequently updated this model to make it specific to dendritic cells (DC) and discovered that RelB is regulated by classical IκBs (IκBα and IκB ) during DC activation and, therefore, integrates the canonical and non-canonical pathways [44]. This computational model is the first to model the kinetics of both RelA and RelB containing dimers.
Despite the successes of these models from the Hoffmann group, it has been argued that they appear to have a number of limitations. The models were calibrated against population-based wet-lab experimentation that used the average of multiple cells rather than single-cell analysis. As such, they did not reflect the large variation in behaviours between individual cells that is inherent to the NF-κB signalling pathway. With this in mind, the work by Nelson et al. [45] built on the early population-based model from the Hoffmann group to incorporate functionality for single-cell dynamics, as seen in their wet-lab results. They used fluorescence imaging of RelA and IκBα to study oscillations in RelA nucleus-cytoplasm localisation in HeLa cells (human cervical carcinoma cells, named after Henrietta Lacks, from whom the cell line was derived [46]) and SK-N-AS cells (human s-type neuroblastoma cells that possess deregulated NF-κB signalling). Their work showed that single-cell time-lapse imaging and computational modelling of RelA localisation showed asynchronous oscillations following cell stimulation that decreased in frequency with increased IκBα transcription. The authors used the ODE-based model developed by Hoffmann et al. and validated against their own single-cell experimental data. Their results demonstrated nuclear-cytoplasmic oscillations of NF-κB during activation and showed that the dynamics of NF-κB translocation are dependent on IκBα transcription, in agreement with work by Carlotti et al. [34], who showed that regulation of NF-κB subcellular localisation is dependent, in part, on the nuclear export function and in part on the cytoplasmic retention function of IκBα.
Lipniacki et al. [47] used the Hoffmann model as a baseline and amended it in three main ways. Firstly, the kinetics of the nuclear and cytoplasmic compartments were refined to take into account differences between their volumes. Secondly, they made use of additional information (i.e., assumptions and considerations from publications by external groups) regarding IκBα interactions within the system; Carlotti et al. [48] and Rice and Ernst [49] suggest that only 10%-15% of the total NF-κB is not complexed to IκBα in a resting cell (the work by Pogson et al. [50] suggesting that 66% of IκBα binds to the cytoskeleton had not taken place yet). Thirdly, they re-estimated the mRNA transcription and translation coefficients through the use of published molecular level data.
Lipniacki et al. utilised a high degree of abstraction and again focused on the NF-κB-IκB signalling module by using 15 ODEs to model the kinetics of IKK, NF-κB, IκBα and the IKK inhibitor A20 within the nuclear and cytoplasmic compartments of cells. They removed IκBβ and IκB and approximated the collective action of all IκB isoforms by IκBα, which is the most active and abundant isoform and also the only isoform that is lethal if knocked out. They also included three new assumptions around the dynamics of the IKK complex: (i) resting cells have a neutral form IKK n ; (ii) stimulated cells have an active form IKK a , which phosphorylates IκBα; and (iii) inhibition of IKK was incorporated into the model, with A20 converting IKK a into the inactive form IKK i , as shown by Krikos et al. [51].
Lipniacki et al., through reference to Yang et al. [52], modified the NF-κB nuclear transport coefficient (from the Hoffmann model) to incorporate a longer time between the IKK-mediated phosphorylation of IκBα within the NF-κB-IκBα complex and the translocation of free NF-κB from the cytoplasm into the nucleus. This modification enabled the model to also take into account the time needed for ubiquitination and proteolysis, within a single coefficient. Additional assumptions were built into this revised model around the IKK complex and concerned: the different binding states of the molecule; the constant synthesis to replace degraded IKK molecules; and the inhibition of IKK by A20, which was also synthesised as part of the model. Their model therefore incorporates two negative feedback loops for NF-κB activation, one involving inhibition of free NF-κB via the binding of IκBα and the other involving the inhibition of IKK by A20. These amendments were made to the existing code (from Hoffmann et al. [35]), and additional scripts to analyse the resulting simulation data were developed in MATLAB. The model was validated against published data of Lee et al. [53] regarding A20-mediated dynamics and the model dynamics of the original Hoffmann model. After parameter fitting, the authors believed that the model successfully reproduced the time behaviour of wild-type and A20-deficient cells.
Following the publication of this augmented computational model, the Hoffmann group (Cheong et al. [54]) reimplemented their model using Cellerator [55], updated their parameter values to take account of the revisions from Lipniacki et al. and extended the model scope to incorporate IKK activation. They complemented this computational work with additional wet-lab experimentation (at the population-level of cells) that focused on constant stimulation with varying TNFα doses. They found that both duration and dose of TNFα have little effect on the duration of the initial NF-κB response and that NF-κB responds sensitively to a wide range of TNFα concentrations. Their in silico experimentation predicted that these signal transduction properties were crucially dependent on the transient nature of IKK activity. Furthermore, they found the dynamics of IKK activity to be non-linear in nature, which they conjecture could be the basis for ensuring robust TNFα-induced NF-κB responses to offset limitations imposed through the diffusion of the ligand within the extracellular environment. More recent work by Tay et al. [56] and Turner et al. [57], who incorporated a number of stochastic processes within their computational models (see semi-stochastic models below), have shown that in contrast to these population-level studies, the NF-κB response is heterogeneous at the single-cell level, with fewer cells responding to lower doses of TNFα, due to the individual cell responses being controlled by a stochastic threshold and, thus, yielding an all-or-nothing response [57].
Finally, the most recent deterministic model that we are aware of was developed by Choudhary et al. [58], who predicted through in silico experimentation (and validated using siRNA knockdown wet-lab experiments) that the canonical and non-canonical pathways are coupled through the action of TNF associated factor 1 (TRAF1) and NF-κB-inducing kinase (NIK). They demonstrated that TNF stimulation (of the canonical pathway) induces TRAF1 expression and that the newly expressed TRAF1 binds to NIK (of the non-canonical pathway) with high avidity. They conjecture that the TRAF1-NIK complex is a central component for cross-talk between the two NF-κB pathways, as the TNF-induced delayed activation of the non-canonical pathway is dependent upon a feed-forward mechanism activated by TRAF1 expression from the canonical pathway.

Semi-Stochastic (Hybrid) Differential Equation Models
Following from these initial ODE models, Lipniacki et al. have published two subsequent increments of their model, by adding stochastic gene activation [59] and then stochastic receptor activation by TNFα [60]; with the wider group also publishing two additional increments, by adding granularity at the single-cell level [56] and incorporating positive feedback through TNFα expression [61]. The first increment continued to use ODEs, but was reimplemented using the MATLAB development tool, focused at the level of a single-cell instead of an averaged population-level, and also introduced a stochastic switch to account for the activity of the genes relating to A20 and IκBα. They did this by using four equations that account for the binding and dissociation probabilities of NF-κB molecules to regulatory sites (in DNA) for A20 and IκBα promoters, the gene transcription into mRNA and the resultant translation into A20 and/or IκBα proteins. This stochastic gene activation facilitates the variability of protein concentrations within simulated cells, which is akin to that seen in biology, e.g., the stochastic switch allows large variances in simulations due to amplifying effects of mRNA synthesis and protein translation. We believe that this model by Lipniacki et al. is the first to have successfully modelled the amplification cascade of the NF-κB signalling pathway following gene activation.
The second increment by Lipniacki et al. used 15 differential equations to model the kinetics of NF-κB signalling, but this time further increased the scope to incorporate the TNFα cell membrane receptor and the enzyme IKK kinase (IKKK), which activates the IKK complex. The TNFα receptor was incorporated into the amplification cascade described above (for the second generation Lipniacki model) through its activation being stochastic. The authors performed manual parameter fitting for this receptor activation, so that 90% of simulated cells are activated in the first 10 min of TNFα stimulation, and claim that predictions from their single-cell level simulation results agree qualitatively with published IKK and NF-κB activity data of the population-level results from Cheong et al. [54]. However, unlike the population-level results of Cheong et al. [54] that suggest that TNFα dose has little effect on the duration of the initial NF-κB response; their in silico experiments show that at a low TNFα dose, only a fraction of cells are activated, but that in these activated cells, the amplification mechanisms assure that the amplitude of the NF-κB nuclear translocation remains above a threshold. Additional in silico experiments showed that low nuclear NF-κB concentration only reduces the probability of gene activation, but does not reduce the gene expression of those responding. They hypothesise that the two effects provide stochastic robustness in responding cells, allowing cells to respond differently to the same stimuli, but causing their individual responses to be unequivocal. This suggests that amplification-saturation dynamics are present within the model, due to the final cell response at the level of NF-κB target genes being approximately equal, regardless of whether a single receptor or 100 receptors were activated. It is currently unknown whether a single activated receptor is sufficient to initiate a response of this magnitude in biological systems, so further wet-lab experimentation is required in parallel to the computational work of this group.
The third increment was developed by Tay et al. [56], who augmented the mathematical model with additional parameters so that it may be used with high-throughput single-cell resolution wet-lab experimental data. The model contains 16 differential equations (some of which are stochastic in nature) and 34 rate constants, 20 of which used fixed constants from published data, and the remaining 14 were manually fitted in order to calibrate simulations against single-cell traces. The most intriguing finding from their wet-lab and in silico experiments, which was also independently discovered by Turner et al. [57] during a similar time period, was that not all cells responded to TNFα and that the fraction of activated cells decreased with decreasing TNFα dose, thus representing the discrete nature of single-cell activation. They also found that early gene expression was not dependent on the intensity of the inducing signal, but instead relied on the high amplitude of the NF-κB, which supported their earlier hypothesis of robust NF-κB responses, even with relatively few cells responding. This model has been used by Fallahi-Sichani et al. [62] as the intracellular basis of a hybrid model, by linking to an existing high-level agent-based model of macrophage responses to Mycobacterium tuberculosis (TB) infection.
They surmise that through manipulating NF-κB-mediated responses (particularly macrophage activation and TNFα expression), you can improve the function of a TB granuloma to contain the infection.
The fourth and final increment from the Lipniacki group was developed by Pekalski et al. [61], who built on the model of Tay et al. [56] by adding the positive feedback associated with the expression of TNFα. Their premise was that the first phase of the innate immune system detects pathogens through membrane and cytoplasmic receptors and that this leads to activation of transcription factors (such as NF-κB) and the production of proinflammatory cytokines, such as TNFα. The secretion of these cytokines then leads to the second phase of the innate immune response in cells that have not yet encountered the pathogen; with these cytokine-activated cells producing and secreting the same cytokine and, thus, propagating the immune response through this positive feedback. They discovered that the introduction of positive feedback changes system dynamics and may lead to long-lasting NF-κB oscillations in wild-type cells and persistent NF-κB activity in A20-deficient cells.
Additional work using hybrid models at the single-cell level has been performed by Ashall et al. [63] and Paszek et al. [41]. Ashall et al. built on the findings of Nelson et al. [45] discussed above. Instead of using the computational modelling in a predictive capacity, Ashall et al. first commenced with wet-lab experimentation and then developed both deterministic and semi-stochastic models to simulate the cellular behaviours. During the wet-lab experimental phase, they exposed individual cells to pulses of TNFα at various time intervals between the pulses, to mimic the pulsatile nature of inflammatory signals. They discovered that lower frequency stimulations generated synchronous translocations of NF-κB across the nuclear membrane of equal magnitude with successive pulses; whereas higher frequency stimulations generated synchronous translocations with reduced magnitude for successive pulses, suggesting that this reflects a failure of the system to reset at higher frequencies. A revised deterministic model, which modified the core network of IKK-NF-κB-IκBα and the A20 inhibitor, which was able to replicate that, was the basis for Lipniacki et al. [47], which was able to replicate the TNFα pulsatile stimulation data at the population-level, using a single parameter set. As the deterministic model was unable to elucidate the heterogeneity of single-cell responses to TNFα pulsatile stimulation, they augmented the model with stochastic processes for three negative feedback mechanisms: transcription of IκBα, transcription of A20 and delayed transcription of IκB . Like the deterministic model, this semi-stochastic model was able to simulate the wet-lab experimental results and predict persistent oscillations after TNFα stimulation, suggesting that stochastic variation due to the delayed transcription of IκB may generate increased cell-to-cell heterogeneity.
Paszek et al. [41] subsequently built on this work of Ashall et al. by augmenting the transduction pathway of the computational model with the kinase IKKK and calibrating the semi-stochastic behaviour to the results of their new single-cell (wet-lab) experimentation into the kinetics of IκBα and IκB activation. Through the use of TNFα and IL-1β cytokines and the synthetic stimulus PMA (phorbol ester differentiation factor phorbol 12-myristate 13-acetate), they discovered a transcriptional delay between IκBα and IκB . Furthermore, they discovered a 45-min delay to transcription of both IκBα and IκB genes and that this delay did not substantially change the average timing or amplitude of NF-κB oscillations at the population-level, but instead affected the single-cell timing of the oscillations to maximise the heterogeneity. They therefore hypothesise that the network topology of the NF-κB signalling system is stimulus-dependent and that the generation of cellular heterogeneity maintains the functional responsiveness of individual cells, to promote tissue robustness.
Through these separate iterations, the models of Lipniacki et al., Tay et al., Ashall et al. and Paszek et al. have now refined the initial Hoffmann model and indeed investigated new areas of the signalling pathway from a computational perspective. We believe that the major achievements are: consideration of the different nuclear and cytoplasmic volumes; introduction of stochasticity regarding cell membrane receptor activation and gene activation, to facilitate modelling at the single-cell level; the use of more accurate data as parameter values; and perhaps, most importantly, the rewriting of the model using MATLAB, which is a more contemporary software development tool for computer scientists.

Agent-Based Models
Although quantitative mathematical models are well established tools for modelling complex biological phenomena, they require an exhaustive set of precise parameters to be specified for each variable. This is fine for small models that align to the dynamics of a few components; however, when the scale is increased to capture a more realistic scope at the system-level [64], they begin to suffer from limitations in accuracy, as data for use in analysing the effectiveness of differential equation models is often unavailable from wet-lab immunological experiments [8]. As such, equation-based methods are unable to fully model system dynamics at the individual component-level and, therefore, suffer considerably from their inability to capture the natural variation inherent to all biological processes. An alternative to these equation-based approaches is to model molecules or cells individually and assign probabilities to each possible interaction or state change through rule-based techniques. The individual component behaviours may then be aggregated up to system-level dynamics, which are then extrapolated in order to make predictions of the system-level behaviours in real biology [65,66].
Interestingly, it would appear that a large amount of data generated experimentally in biology actually accumulates in an object-oriented manner. The reductionist approach, which endeavours to look at systems using the smallest indivisible unit, is analogous to looking for "objects" within nature. Object-oriented approaches to modelling, therefore, provide a useful formalism for constructing computational models by designing systems from a bottom-up perspective and organising information around individual objects [67]. Agent-based modelling and simulation (ABMS) builds on the object-oriented paradigm, with the key enhancement being that an agent is active rather than passive and that ABMS has multiple threads of control. Macal and North [68] nicely describe this by stating that the "fundamental feature of an agent is the capacity of the component to make independent decisions". The agent-based approach therefore allows systems to be modelled at the resolution of individual components and allows us to model at the level of the individual molecule the spatio-temporal aspects of system dynamics, such as translocation of components between cell compartments, and also simulate their individual position in three-dimensional space.
Following from the early successes with equation-based computational models, Pogson et al. [69] took the modelling of the pathway's stochastic behaviour one step further, by developing an agent-based model of the NF-κB signalling module. This computational model was developed in MATLAB and utilised the concept of communicating X-machines [70] to represent the individual agents and their associated interactions. They chose to break away from the pattern of using differential equations, which relied on the assumption that a cell is homogenous with respect to its chemical constituents. We agree with this stance and believe that such an assumption will not be valid, due to the cells' internal compartmentation and the non-uniform distribution of key molecules [9]. An important aspect of communicating X-machines is that each agent has memory, which in this instance, holds the current physical location and current state and may also contain a set of randomised parameter values (from a given distribution of suitable parameter values), to further instil stochasticity into the model. Pogson et al. used a set of autonomous agents within their model that correspond to: NF-κB dimer, IκBα, cell membrane receptor, nuclear importing receptor and nuclear import-export receptor. Furthermore, they used a greater level of detail regarding translocation of molecules across the nuclear membrane and added a three-dimensional spatial dimension using continuous space to the model, which the differential equations used by Hoffmann et al. and Lipniacki et al. are unable to represent, and to a higher degree of granularity than the partial differential equations (PDEs) used by Terry and Chaplain [71] and Ohshima et al. [72] could represent.
The mathematical underpinning of this agent-based model is based on generalised biochemical reaction kinetics of two substrate molecules interacting and forming a product molecule. Graphical visualisation of the reactions allows us to view the dynamics of the system over time, for example the increase of product molecules (e.g., NF-κB-IκBα complex) and the associated decrease in substrate molecules (e.g., free NF-κB and free IκBα) during an inhibitory reaction. As the relevant components are located within three-dimensional space (X, Y, Z co-ordinates), the authors built in an interaction radius to each molecule to ensure that a chemical reaction only occurs when the requisite substrate molecules come within a certain distance of each other. They abstracted the interaction radius of molecules to be a sphere; clearly, a number of biochemical components have orientations, through, for example, polarity and non-symmetrical quarternary structure; however, we believe this to be an appropriate approximation to biology for the purposes of the model.
Simulations (using the model) were run as a discrete-event system, with each iteration of the simulation run representing a defined period within time, i.e., constant time step. The authors validated the model against single-cell analysis wet-lab data of Carlotti et al. [34,48] and Yang et al. [73], who used fluorescent tagging and confocal microscopy to observe IκBα degradation and NF-κB shuttling dynamics. Successive time-step iterations were therefore able to relate to the temporal observations of single-cell analysis experiments, and the model therefore reflects the discrete stochastic nature of sub-cellular events.
Pogson et al. later updated the original model and performed in silico experimentation to predict internal cell structural components in regulating the NF-κB signalling pathway [50]. This revised model made the IKK complex explicit and introduced the transcription and translation of the three IκB isoforms (IκBα, IκBβ and IκB ), as per the model of Hoffmann et al. IκBβ and IκB were included in the model to introduce more accurate stochasticity in the inhibition of NF-κB and reduce the overall probability for NF-κB-IκBα complex formation. As with the 2006 model, this augmented model was also validated against single-cell analysis data [34,73] and sensitivity analysis demonstrated a narrow range of acceptance for IKK levels, robustness to NF-κB and IκBα levels and importantly greater sensitivity to IκBα than for NF-κB. Furthermore, simulation results yielded a negative feedback loop consistent with the models of Hoffmann et al. and Lipniacki et al.
Through in silico experimentation, the authors demonstrated an increase in NF-κB nuclear translocation at a 1:1 NF-κB:IκB ratio (which is to be expected) and focused on the cytoskeletal structural component of the cell to which IκB can bind. Simulations of cytoskeleton-IκB interactions yielded a maximal NF-κB-IκB complex formation at a 1:3 ratio NF-κB:IκBα. Without actin, the maximal was reached at a 1:1 concentration ratio. The authors hypothesise that this key role for cytoskeleton-associated interactions sustains optimal pathway regulation by adjusting the NF-κB-IκB complex formation at the steady-state and controlling negative feedback following activation of the signalling pathway. Subsequent wet-lab experiments (performed by the authors) demonstrated a ratio of actin-bound to free IκBα of 2:1 in unstimulated cells, indicating that two-thirds of IκBα may be bound to the cytoskeleton at steady-state.

Peer-Validation of These Computational Models
Adverse effects of over-transfection were initially observed by Carlotti et al. [48] and Yang et al. [73], who compared the control of exogenous, tagged intermediates with results from immunocytochemical experiments to identify transfection limits in regards to cell-to-cell variations and regulation of NF-κB. This risk of over-transfection was initially leveraged by Barken  We also believe that continued refinement and augmentation of these existing models will in all likelihood yield a computational model that will be able to accurately predict observed biological phenomena across both single-cell and averaged population levels, and we suggest that the recent work of Tay et al. [56] and Turner et al. [57] has made an important contribution towards this end. They have both shown through single-cell analysis with TNFα and associated stochastic modelling that NF-κB activation is heterogeneous and is a digital process, with fewer cells responding at low doses. Additionally, Turner et al. also found that cells display analogue processing to modulate system dynamics, by controlling NF-κB peak intensity, response time and the number of oscillations.
Hayot and Jayaprakash [78] incorporated the Gillespie algorithm [79] into the model for the study of cell-to-cell variability, along with the replacement of the quadratic (second order) nuclear NF-κB-induced synthesis of IκBα with a linear (first order) term, following the reasoning of Nelson et al. [45] and Lipniacki et al. [47]. This augmented model found the fluctuations of NF-κB to agree with Hoffmann et al. [35] and predicted that these fluctuations are small when transcription is strong (as used within the Hoffmann model), whether transcription is linear or quadratic, but that intrinsic fluctuations can be strong and, indeed, as significant as any extrinsic noise if the rate of promoter binding is low.
Other groups have used models to validate their in-house wet-lab experimental data, including: Ihekwaba et al. [80,81], who performed a parameter refit for IκBβ and IκB , discovered, through using sensitivity analysis, that the model is particularly sensitive to parameters relating to IKK and IκBα and that overexpression of IκBα in transfected cells would yield a minimal perturbation to the system; Mathes et al. [82], who added IκBα degradation functionality for both IKK-mediated degradation of NF-κB-bound IκBα and IKK-independent degradation of free IκBα; and Wang et al. [83], who adapted the model to simulate single-cell data using short and strong pulses of TNFα and found that weak pulses provide a rich variety of non-linear temporal signalling dynamics, which may account for the diversity of expression patterns for the multitude of target genes in the signalling pathway.
Similarly, Yue et al. [84], who performed a sensitivity analysis, focusing in particular on local sensitivities through bifurcation/pairwise studies; Joo et al. [85], who performed global sensitivity analyses using orthogonal array sampling, Latin hypercube and k-means clustering and confirmed that the IKK-NF-κB-IκBα-A20 interaction network is the key portion of the wider pathway (as suggested by Lipniacki et al. [47,59]), as these components have critical kinetic rate variables for IKK activation; and Ashall et al. [63], who performed wet-lab experimentation by treating cells with short pulses of TNFα to mimic pulsatile inflammatory signals and, through in silico experimentation on a model, augmented with various inhibitors, predicted that negative feedback loops regulate both the resetting of the system and cellular heterogeneity, through limiting the reactivation of IKK.

Minimal Models
The models discussed above have represented mathematical and computational abstractions of the major components of the NF-κB signalling pathway. It has recently been argued by a number of different groups, however, that the complex nature of these models makes them computationally expensive to run and also requires complicated analysis techniques in order to make inferences from simulation results. As such, these groups have developed what have been termed minimal models, which are able to replicate the majority of the phenomenological behaviours of the more complex computational models, but which use the minimum number of equations possible. We therefore discuss within this section the four minimal models from Krishna  The first minimal model after that of Carlotti et al. [34] was developed by Krishna et al. [86] and focused on a small core network of the pathway that drove oscillatory behaviour. Through reducing the model of Hoffmann et al. down to a core feedback loop of three coupled ODEs, they were able to ascertain the minimal model required to generate oscillations. This simplified model was validated against the work of Hoffmann et al. and Nelson et al. and was able to simulate: the sustained oscillations obtained with only the IκBα isoform of IκB; the damped oscillations in wild-type cells that included the other forms of IκB (e.g., IκBβ and IκB ); the spikiness of nuclear NF-κB and the asymmetry of cytoplasmic IκB oscillations; and the phase difference between NF-κB and IκB. Their key finding was that saturation degradation of cytoplasmic IκB in the presence of IKK was crucial for oscillatory behaviour, because it sets an upper limit within the system for the degradation rate and, thus, allows IκB to accumulate and remain in the system for longer than with linear degradation rates. They conjecture that this effectively introduces a time delay into the negative feedback loop, which is known to generate oscillations.
The authors were subsequently involved in the development of another minimal model (Yde et al. [87]). Here, they modelled the amplification of the immune response, following cytokine-mediated activation of the NF-κB pathway, and the positive feedback that occurs through upregulation of cytokine expression (as previously modelled by Werner et al. [88]). The ODE model was developed at the tissue-level (cell population level), with the detailed network dynamics of the NF-κB signalling pathway within individual cells abstracted away. In fact, they discovered that tissue-level immune responses were able to emerge within the minimal model through using only three (high-level) variables, relating to: NF-κB, a generic regulator and a generic inhibitor. The model predicts that cytokines produced by the stimulated cell(s) at the site of infection diffuse away from this primary infection site and trigger the transient response of the NF-κB pathway for the production of cytokines in neighbouring cells. They conjecture that this generates a propagating wave of NF-κB induction and cytokine production throughout the infected tissue.
Like the minimal model of Krishna et al. [86] above, Longo et al. [89] also used the Hoffmann model as the basis for developing a minimal model of negative feedback within the NF-κB signalling pathway. This new model used a single delayed compound reaction to replace the cascading reactions that had been utilised previously [35]. They used this to investigate the underlying mechanisms involved with oscillations and the dampening of oscillations that emerged through the two negative feedback loops of the inhibitors IκBα and IκB . Through investigating the oscillatory responses using the single negative feedback associated with IκBα, they discovered that both the frequency and decay rate of the oscillations are highly dependent on the internal parameters of the core network of the signalling module, but are not sensitive to extracellular stimuli levels. They suggest that the oscillatory frequency within the system can therefore not be encoding information about the stimulus and conjecture that stimulus-specific gene expression is therefore unlikely to be determined by the stimulus-specific frequencies of NF-κB oscillations; and therefore, may involve amplitude modulation over time. Furthermore, upon introduction of IκB into the model, the authors were able to reproduce the findings of their previous work [40], where the oscillations caused by persistent stimulation were now significantly dampened. They suggest that the second negative feedback mechanism (IκB ) may have evolved to produce dampening of the oscillatory behaviour of the first feedback mechanism (IκBα).
An alternative approach has been taken by Zambrano et al. [90], however, whereby the ODE-based model is developed according to the different layers of processes within the cell environment. They have used a three-layer approach, with the first layer representing the transcription and translation of new IκBα molecules; the second layer representing the NF-κB signalling module and its negative feedback by IκBα; and the third layer representing the commencement and propagation of signal following the recognition of extracellular stimuli at the cell membrane, down to the activation of IKK. Their model was able to replicate both the heterogeneity of the system seen in the models of Tay et al. [56], Paszek et al. [41] and Sung et al. [76], along with the spiky oscillations shown by Krishna et al. [86]. As such, the authors believe that the community should make more use of minimal models as a basis for investigating the underlying mechanisms of the NF-κB signalling pathway and that the reliance on complicated computational models should be minimised, because they incur the risk of being overfitted to the specific wet-lab experiments from individual groups and, thus, may not be used by the wider community, who use different in vitro models.
Finally, along with the explicit development of minimal models by individual groups, recent work by West et al. [91] has focused on the development of algorithmic methods to enable a principled approach for reducing existing differential equation-based models into new minimal models. Their algorithmic approach is based on quasi-steady-state approximations and produces a set of ranked variables according to how quickly they approach their momentary steady-state. Once ranked, the developer may then eliminate variables at each step within the related equation-based model, whilst preserving the system-wide dynamics. The authors tested their reduction algorithm on the two feedback (IκBα and A20) models of Ashall et al. [63] and the original ODE model by Krishna et al. [86] (upon which the minimal model was extracted) and report that the system dynamics that emerge from the reduced models compares favourably with the original, more complex models.

Discussion and Perspectives
The complexity of the innate immune system, and in particular the inflammatory process, has been difficult to fully investigate using reductionist and linear approaches alone, since it is characterised by non-linear kinetics, as well as numerous feedback loops. Several signalling pathways relating to the combating of various pathogenic infections converge on NF-κB activation, resulting in a highly complex regulatory system for the innate immune response. Since the discovery of NF-κB in 1986 [12], there have been over 49,000 journal articles published (figures from PubMed on March 5, 2014). We believe the NF-κB signalling pathway to be a good candidate for research that follows a systems biology approach to further our understanding. Computational modelling to date has successfully complemented traditional wet-lab techniques to generate hypotheses on the mechanistic behaviour of the various components of the signalling pathway (see Table 1 for the major computational modelling advancements). As such, these in silico-derived predictions are helping to further our understanding of in vitro system dynamics, in particular how signals received by the IKK complex are propagated down the signal transduction pathway through activation of the NF-κB signalling module. This computational experimentation has given rise to the belief that IKK regulatory mechanisms may represent sensitive clinical targets in diseases with aberrant innate immune system activity. There appears to be an opportunity to expand the scope of these current models, however, to also encompass upstream events prior to IKK activation, thus incorporating a greater degree of granularity for cell membrane receptor complexes, and to expand the scope to additional gene transcription products. This would allow us to simulate the induction of the pathway through various extracellular signals and to capture activation of transcription for target genes (other than IκB isoforms) that are relevant to the inflammatory response, for example various cytokines. Table 1. The top 20 computational models (in chronological order) that we believe have facilitated increased understanding of the NF-κB signalling pathway. Along with the year of journal article publication, we have also documented: the modelling paradigm used (deterministic, semi-stochastic, or agent-based), the specific extracellular stimulus (if any), the specific cell type (if any), the hierarchical level of the model (single-cell or population), the pathway components that were explicitly modelled, and the key advances of the model.  NF-κB has wide-ranging effects controlled by a complex regulatory network of inhibitors and co-activators. Understanding the mechanisms that control NF-κB activation/cellular signalling is important for exploiting therapeutic approaches to treat human disorders due to its dysregulation. Specific targets for therapeutic agents could be the transcription factor itself or any of the associated components within the pathway, including protein kinases and the IκB inhibitors themselves. Gilmore and Herscovitch [92] reviewed the known inhibitors of NF-κB, which may provide a basis for future research regarding computational modelling and subsequent pharmacological intervention.
In conclusion, we believe that computational models tightly coupled to wet-lab experimental analysis will be indispensable to furthering our understanding of the NF-κB signalling pathway. Furthermore, with the advances in software modelling technologies and computational power over the past decade, we believe that the timing is right for the development of agent-based models of the NF-κB signalling pathway on a scale that has not been seen before. New technologies, such as Java Mason [93], which allows large agent-based simulations to be developed, which can be manipulated in real-time during simulations, or FLAME (Flexible Large-Scale Agent-Based Modelling Environment) [94,95], which has recently been used to model the European economy, and the GPU version of FLAME [96], allow massively parallel agent-based simulations to be run. Along with these advances in technology, recent developments in cell biology have also advanced our corresponding understanding the pathway. We believe a number of key questions may be answered through the development of large-scale computational models of NF-κB. From a cell biology perspective, these are: (i) what can in silico experimentation tell us about the relative roles of the intermediate components within the signal transduction events; (ii) what can in silico experimentation tell us about the various receptors, co-receptors and adaptor proteins and their role(s) in signal transduction events; (iii) what can in silico experimentation tell us about the dysregulation that can occur in diseased states, and how can we perturb the system back to a healthy state; and (iv) what can in silico experimentation tell us about the cross-talk that occurs when various extracellular signals converge on the NF-κB signalling module. There are also a number of key questions from a modelling perspective: (i) does the abstraction level affect the accuracy of simulation-driven predictions; (ii) does the resolution level (e.g., number of agents) affect the accuracy of simulation-driven predictions; (iii) does the simulation platform (different technology, e.g., equation-based versus agent-based) affect the accuracy of predictions; (iv) what are the relative merits of averaged population data versus single-cell data for the calibration and validation of computational models; (v) what are the advantages and limitations of using massively parallel computing architectures, for reproducing the large-scale variation (as seen in biological systems) into computational simulations; and (vi) what can inference techniques (e.g., Bayesian inference) tell us about the relative merits of the existing computational models and their associated parameter value sets?