 Research
 Open access
 Published:
Propagation of kinetic uncertainties through a canonical topology of the TLR4 signaling network in different regions of biochemical reaction space
Theoretical Biology and Medical Modelling volumeĀ 7, ArticleĀ number:Ā 7 (2010)
Abstract
Background
Signal transduction networks represent the information processing systems that dictate which dynamical regimes of biochemical activity can be accessible to a cell under certain circumstances. One of the major concerns in molecular systems biology is centered on the elucidation of the robustness properties and information processing capabilities of signal transduction networks. Achieving this goal requires the establishment of causal relations between the design principle of biochemical reaction systems and their emergent dynamical behaviors.
Methods
In this study, efforts were focused in the construction of a relatively well informed, deterministic, nonlinear dynamic model, accounting for reaction mechanisms grounded on standard mass action and Hill saturation kinetics, of the canonical reaction topology underlying Tolllike receptor 4 (TLR4)mediated signaling events. This signaling mechanism has been shown to be deployed in macrophages during a relatively short time window in response to lypopolysaccharyde (LPS) stimulation, which leads to a rapidly mounted innate immune response. An extensive computational exploration of the biochemical reaction space inhabited by this signal transduction network was performed via local and global perturbation strategies. Importantly, a broad spectrum of biologically plausible dynamical regimes accessible to the network in widely scattered regions of parameter space was reconstructed computationally. Additionally, experimentally reported transcriptional readouts of target proinflammatory genes, which are actively modulated by the network in response to LPS stimulation, were also simulated. This was done with the main goal of carrying out an unbiased statistical assessment of the intrinsic robustness properties of this canonical reaction topology.
Results
Our simulation results provide convincing numerical evidence supporting the idea that a canonical reaction mechanism of the TLR4 signaling network is capable of performing information processing in a robust manner, a functional property that is independent of the signaling task required to be executed. Nevertheless, it was found that the robust performance of the network is not solely determined by its design principle (topology), but this may be heavily dependent on the network's current position in biochemical reaction space. Ultimately, our results enabled us the identification of key rate limiting steps which most effectively control the performance of the system under diverse dynamical regimes.
Conclusions
Overall, our in silico study suggests that biologically relevant and nonintuitive aspects on the general behavior of a complex biomolecular network can be elucidated only when taking into account a wide spectrum of dynamical regimes attainable by the system. Most importantly, this strategy provides the means for a suitable assessment of the inherent variational constraints imposed by the structure of the system when systematically probing its parameter space.
Background
Normal and abnormal cellular states represent macroscopic behaviors emerging from intricate dynamical patterns (either transient or stationary) of biochemical activity. These are sustained by a complex web of reaction mechanisms that play the role of information processing systems, generically referred to as signal transduction networks [1ā3]. In other words, these networks represent the dynamical systems that instruct cells to enter into specific regimes of biochemical activity, which ultimately determine the universe of functional states accessible to the cell, such as differentiation, apoptosis, cell division, etc. [1ā3]. Operatively, functional regimes of biochemical activity within a cell are basically accomplished via direct proteinprotein interactions and enzymecatalyzed reactions (i.e. phosphorylation, RNA synthesis, etc.) triggered in response to either internal or external stimuli [3, 4].
The spectrum of functionalities that a signal transduction network can potentially perform is inherently constrained by its design principle [5, 6], which encapsulates a series of aggregated components involving diverse regulatory schemes and biochemical reaction rules modulated quantitatively via internal reaction parameters. This structurefunction puzzle has motivated considerable research efforts in the last decade aimed at elucidating possible mechanistic bases of fundamental emergent properties such as robustness, evolvability and epistasis, of highlymodular regulatory systems [7ā13]. Importantly, the investigation of the robustness properties of a signal transduction network requires heavy emphasis to be made on two fundamental aspects of the underlying reaction mechanism: an observable/quantifiable dynamical feature (either transient or stationary) of the system, and one or several perturbable parameters directly or indirectly involved in the development of the system's feature being studied. For instance, important quantitative dynamical features of signal transduction networks have been proposed as suitable targets for assessing their robustness properties in the face of random changes in internal reaction parameters [14, 15]. Sources of perturbations impinging upon such parameters may stem from environmental vicissitudes (temperature, pH, etc.), genotypic variation or intrinsic fluctuations (molecular noise) [16, 17].
Recently, several computational studies have yielded interesting numerical evidence supporting the idea that the robustness properties of highlydimensional biochemical reaction networks may be strongly dependent on three fundamental aspects: i) the reaction topology (network architecture) [7ā9], ii) the system's current position in parameter space [18ā20], and iii) the dynamic nature of the trajectories displayed by the reaction species involved [13, 20ā22]. The robustness properties of a biomolecular network are typically assessed by means of standard sensitivity analysisbased approaches implementing both local and global perturbation methods [18, 23ā27]. Robustness is usually assessed with respect to either observable or hypothetical stationary states and transient dynamics of just few reaction species in the network [24, 28, 29]. However, a complementary quantitative approach to studying the robustness properties, as well as information processing capabilities, of a complex reaction network should provide the means for assessing the extent to which the full dynamical behavior of the system is reproducible under, for example, kinetic uncertainties. This is because a reaction network may be coupled dynamically in unexpected ways to other important subsystems not included in the model [11, 30], whereby biochemical information exchange among cellular processes can take place in parallel. Under these considerations, we thus believe that general properties of a canonical biomolecular network could be revealed under the following methodological strategies. Firstly, a large ensemble of disparate, but biologically plausible dynamical trajectories attainable by the network should be tested for general robustness properties in the face of random perturbations impinging upon the whole set of reaction parameters; that is to say, the overall robust performance of the network should be evaluated in widely scattered regions of its accessible parameter space. Secondly, the reproducibility of particular ouputs (i.e. experimentally reported wildtype transcriptional readouts) should be assessed in different regions of the accessible parameter space via both local and global perturbation strategies. Addressing these points would pave the way to gaining general insight into systemslevel features of the complex reaction mechanisms endowing the cells with the potential to reach a wide spectrum of robust behaviors.
In this study, efforts were focused on a comprehensive and unbiased statistical assessment of the robustness properties and information processing capabilities of a canonical reaction topology underlying TLR4mediated signaling events. This signaling network is temporally deployed in inflammatory cells (i.e. macrophages) in response to external stimuli. We constructed a deterministic, nonlinear dynamic model of this reaction topology, using an informational basis retrieved from a series of previous computational studies and review papers providing important clues about mechanistic reaction steps involved in the process (see the Results and Discussion section below). We adopted this signaling network as our model system mainly because this functional module plays a crucial role in the development of innate immune cellular responses ([31ā37]). For instance, Tolllike receptors recognize conserved pathogenassociated molecular patterns such as lipopolysaccharide (LPS), which results in the triggering of both microbial clearance and the induction of immunoregulatory chemokines and cytokines. Here, we centered our attention specifically on the immediate cellular response, in macrophages, triggered by the rapid activation of the canonical MyD88dependent and TRIFdependent reaction cascades upon LPS binding to TLR4. We probed the robustness properties and information processing capabilities of this canonical network in different points distributed across diverse regions of the biochemical reaction space. Importantly, the behavior of the network in a given region of the biochemical reaction space was selected so that it was congruent with a hypothetical, but biologically plausible dynamical regime of molecular activity (see below). Global (nonorthogonal) and local (orthogonal) perturbation strategies were implemented as a means of systematically exploring the biochemical reaction space inhabited by the network. Critically, reaction parameters were subjected to random perturbations without a priori knowledge on their relative importance for the network in the accomplishment of a given signaling task. Our extensive numerical analyses permitted us the identification of global and particular variational constraints in the network. This was achieved by means of a detailed characterization of some statistical regularities on the dynamical performance of the system under kinetic uncertainties (i.e. random fluctuations in internal reaction parameters). Overall, our simulation results provide convincing numerical evidence supporting the following idea: a canonical reaction mechanism underlying TLR4mediated signaling events is endowed with the intrinsic capacity to perform information processing in a robust manner, which is remarkably independent of the signaling task required to be executed. Nevertheless, our statistical analysis indicate that the robust performance of the network is not solely determined by its architecture (topology), but this may be strongly conditioned by the network's current position in biochemical reaction space. Ultimately, our simulation results provide interesting mechanistic insigths into structurefunction relationships in the TLR4 signal transduction network, which enabled the identification of plausible rate limiting steps that most effectively control the performance of the system under diverse dynamical regimes.
Information processing and biochemical reaction space of the signal transduction network
To avoid any confusion or controversy regarding well stated systems biology concepts on cell signaling processes, it is important to make clear our notion of a signal transduction network as an information processing system, mainly because this may differ considerably from previous conceptualizations. Nevertheless, we believe our conceptualization provides a complementary view of the issue. For example, the notion of information processing applied in the context of intracellular signaling has traditionally been limited to the mechanistic explanation of how cellular behaviors are induced via the decodification, and subsequent intracellular propagation, of time variant/invariant physicochemical signals provided by extracellular stimuli (see for example [6, 38ā43]). Our intent here was to extend the scope of this notion, making it more suitable for systemslevel robustness analysis of signal transduction networks. Our rationale focuses on the following arguments. Given that the emergence of central cellular behaviors relies heavily on the robust performance of signal transduction networks, it follows that the information processing capabilities of these systems are primarily dependent on internal reaction parameters. In general, such parameters exhibit a natural tendency to behave like a set of random variables, resulting mainly from thermal fluctuations in the cell environment, and mutational perturbations in the genetic encoding of the system. Arguably, the internal reaction parameters of a signaling network stand for repositories of kinetic information that collectively define a biochemical reaction space inhabited by the system. Such a reaction space becomes an essential source of information carefully coupled to extrinsic stimuli that turn out to be processed according to the set of reaction rules encoded in the architecture of a signal transduction system, from which a proper cellular phenotype (i.e dynamic protein activation profiles and/or gene expression patterns) is calculated (see Figure 1). Ideally, these should represent the basic tasks any information processing system, such as a signal transduction network, is expected to accomplish in a robust fashion. Under these considerations, it should be clear that we equate robust information processing capabilities of a signaling network with its capacity to reproduce particular (reference) dynamical trajectories of biochemical activity under random perturbations in its internal reaction parameters. Importantly, this is assessed here via standard metrics aimed at evaluating discrepancies between dynamical trajectories, and by means of rigorous statistical analysis (see the "Models and computational framework" section below). Our methodology can thus be seen as a coarsegrained strategy to assessing the information processing capabilitites of a complex reaction network, when monitoring the propagation of kinetic uncertainties throughout the system. This represents an alternative framework to that recently proposed methodology relying on Shannon's entropy (see [44]). Interestingly, that framework conceives a signaling network as a "communication channel", for which the associations between inputs and outputs result from a decomposition of their mutual information into different components.
Methods
Canonical reaction topology underlying TLR4mediated signal transduction events
Within a rather short time window, LPS binding to TLR4 triggers two major intracellular signaling events rapidly propagated through the MYD88dependent and TRAMdependent reaction cascades, which display extensive crosstalking (see Figure 2). Activation of the MYD88dependent cascade leads to induction of proinflammatory cytokines such as TNFĪ± by means of JNK, p38, NFĪŗ B and ERK; whereas the TRAMdependent cascade predominantly induces the expression of chemokines such as the IP10 protein encoded in the Cxcl10 gene, via the interferon regulatory factor (IRF) [45]. A relatively limited number of existing dynamic modeling studies focus specifically on TLR4mediated signal transduction. For example, pioneering simulation works have provided interesting mechanistic insights on diverse kinetic phenomena observed during temporal deployment of this signal transduction network, such as time delay responses [46], signaling flux redistribution [47], and preconditioning behavior [48, 49]. Based upon the information provided by these theoretical studies and the data reported in recent review articles about key architectural features of this signaling network (see for example [31ā37]), we assembled a wellinformed mathematical representation of the complex web of biochemical reactions that are likely to sustain the information processing capabilities of this signal transduction system. Our modeling framework is grounded on ordinary differential equations incorporating first and second order reactions for representing intracellular signaling fluxes, as well as Hilllike saturation kinetics accounting for highly nonlinear reaction schemes taking place at the level of ligandreceptor interactions and transcriptional activation (see "Models and computational framework" section below, and Additional file 1 for a detailed description of the mathematical structure of the network model). The total number of reaction species modeled amounts to 76, including a TLR4 in both a susceptible and an activated form, MYD88 and TRAM adapters along with their associated molecules, hypothetical intermediates upstream to TRAM which have been inferred computationally in [46, 47], intermediate and effector kinases (i.e. MKK4/7, JNK, MKK3/6, p38, TpL2, MKK1/2, ERK), the associated and dissociated forms of NFĪŗ B and IĪŗ B, and two important mRNAs transcribed from the TnfĪ± and Cxcl10 proinflammatory genes (see Figure 2). We also assumed a time variant concentration of LPS following an exponential decay profile as an alternative hypothesis to that simulated intrinsically stable dynamic regime of LPS proposed in a recent study of TLR4 activation kinetics ([48]). Nuclear export and import dynamics from the cytoplasm of some reaction species were modeled via simple first order kinetics, hence, volumedependent scaled coefficients of transport were neglected for simplicity purposes. Moreover, within the narrow time window simulated, our modeling framework assumes that simple first order reaction kinetics govern dephosphorylation processes. In this way, dephosphorylation of a substrate was only dependent on its own concentration and the dephosphorylation rate. Furthermore, we lumped together into single reaction steps multisite phosphorylation processes, which might not represent key rate limiting steps in the cascades included in our model. We therefore have equated multisite phosphorylation steps with full kinase activation, which might constitute a truly rate limiting step during signal processing. It is also worth saying that an explicit mathematical representation of the dynamics of ATP was not considered; instead, we assumed it to be in a steady state. This is standard practice in kinetic modeling and is implemented for simplicity purposes. Our mathematical representation of the whole reaction scheme defines a multidimensional biochemical reaction space encompassing 116 kinetic coefficients (axes), including transition rates between receptor states (susceptible ā activated), production and degradation rates of receptors, association/dissociation rates among intracelular molecular species, phosphorylation/dephosphorylation rates, nuclear import/export rates, maximal transcriptional rates, transcriptional efficiencies, MichaelesMenten constants, cooperative coefficients, and mRNA degradation rates. Reaction kinetic values for this signaling system have so far proven extremely difficult to assess under well controlled experimental conditions. Therefore, our massive amounts of computationally predicted values of internal reaction parameters for this signaling network might provide a glimpse on the kinetics of the system under different cellular states. Moreover, despite obvious simplifying assumptions about the intricacies of the reaction steps involved, our mathematical representation captures core design principles of the signal transduction network. This is because our model was validated with dynamic experimental data (time courses) from wildtype target transcriptional readouts, which have been shown to be actively modulated, in quantitative terms, by the reaction cascades accounted for in our proposed scheme (see below). Critically, our simulated time window was limited to an interval spanning 120 minutes, a time scale during which critical transient transcriptional readouts are realized as a result of rapidly mounted innate immune responses ([47]). Furthermore, the transient features exhibited by the network during such time period emerge primarily as a consequence of intrinsic processes guided by the intracellular regulation of TLR4 signaling in response to LPS. This is opposed to those extrinsic processes triggered by autocrine and paracrine stimuli provided by antiinflammatory cytokines (i.e. IL10 and TGFbeta), which entails the temporal deployment of complex regulatory schemes such as (+/) feedback control. Presumably, within the narrow temporal window of TLR4 activation in response to LPS stimulation, which is the focus of our modeling framework, the initial signaling phase might not be heavily dependent on the complex feedback dynamics that are subsequently displayed by the NFĪŗ B regulatory module [50]. Such dynamics, instead, should play a major role in reliable control of a delayed (secondary) signaling phase in response to LPS stimulation (see for example [51]). Interestingly, the presence of two signaling phases in this crucial immune cellular process might represent very distinct episodes of signaling fluxes, carrying particular information, that differentially modulate in quantitative terms the transcriptional readout of specific gene batteries.
Results
General robustness properties of the signal transduction network in different regions of the biochemical reaction space
Our first round of numerical experiments was designed with the main goal of exploring the intrinsic robustness properties of the whole integrated reaction network. We computationally reconstructed a rather limited ensemble of 100 different signaling regimes or dynamical trajectories (i.e. the set of 76 individual temporal profiles for the reaction species modeled, which is associated with a given point in parameter space) attainable by the network (see Figure 3). We randomly explored the parameter space looking for solutions in which some reaction species undergoing, for example, covalent modifications (i.e. phospho/dephosphorylation) displayed particular dynamic features similar to previously simulated, and experimentally reported, signaling outputs. Specifically, we focused on trajectories displaying biologically plausible dynamical signatures, such as sustained and transient dynamics of molecular activity with identifiable signaling peaks in some cases. Our simulated reference trajectories were thus required to match, at least qualitatively, distinct signaling outputs previously reconstructed computationally from experimental data (see for example [14, 15, 28]). Under these considerations, such an ensemble of reference trajectories can be thought of as being congruent with a plausible spectrum of cellular states attainable by, for example, a macrophage, which may be a natural operative condition (i.e. phenotypic plasticity) of many types of immune cell lineages ([52]). Alternatively, such an ensemble of dynamical trajectories can be seen as a set of widely scattered points in the multidimensional biochemical reaction space (see Figure 4), with some points being closely related and defining small neighborhoods in biochemical reaction space. As noted above, we randomly explored the parameter space according to a previously defined range of variation assigned to each reaction parameter (see Additional file 1 for a detailed description of parameter ranges); ranges of variation were constrained based on previous simulation results obtained from random scrutiny of the parameter space (personal observations, data not shown), and biological intuition. Moreover, each reference dynamical trajectory was propagated from a particular set of intital conditions (see Additional file 1 for a detailed description), which were also constrained based on previous simulation results (personal observations, data not shown) and biological intuition. Initially, thousands of simulated trajectories were carefully monitored both manually and systematically in order to assemble our final ensemble of biologically plausible dynamical trajectories. Comprehensive statitistical analyses were performed over our limited ensemble of reference trajectories. Ultimately, by following this computational methodology, we were able to conduct a series of well controlled in silico experiments that allowed us to probe the intrinsic robustness properties of the network, under different hypothetical scenarios of biochemical activity. We implemented a global (nonorthogonal) perturbation scheme, also known as multiparametric sensitivity analysis (MPSA) (see the "Models and computational framework" section below). This computational methodology provides the means for conducting efficiently systematic rounds of perturbations in each of the 100 reference points (reference parameter configurations) distributed throughout parameter space. Each reference parameter configuration was subjected to a round of 5000 simulataneous perturbations; that is to say, 5000 newly assembled parameter configurations surrounding each reference point in parameter space were generated. To do this, we first performed uncertainty analysis consisting of Monte Carlo simulations based on the efficient Latin Hypercube Sampling (LHS) scheme, followed by sensitivity analysis, which allowed the identification of those reaction parameters most critically involved in the global performance of the reaction network (see [23, 53, 54], and the "Models and computational framework" section below). Importantly, under this framework the robust information processing capabilities of our model reaction network were properly evaluated by means of a detailed statistical analysis of the system's global sentivities. We analyzed the distributions of the D statistics calculated from KolmogorovSmirnov (KS) tests performed by means of the afortmentioned perturbation approach. Briefly, a KS test is intended to evaluate the global sensitivity of the system's output with respect to perturbations targeting individual parameters. This test specifically provides the means for evaluating the cumulative frequency of the observations (parameter values) as a function of class, and for calculating the maximum vertical distance between cumulative frequency distribution curves for m acceptable and n unacceptable cases of any given parameter Īø_{ j }(see the "Models and computational framework" section below). Figure 5 illustrates a series of box plots summarizing the overall statistical tendency of the D values calculated for each reaction parameter of the network model, over the ensemble of 100 dynamical trajectories that were systematically perturbed. Here, it is worth noting that for each perturbation study, the perturbed signaling trajectories were compared only with a corresponding reference trajectory; being such a trajectory a member of the ensemble of 100 trajectories analyzed. In general, our analysis indicates that the network is capable of reproducing reference dynamical trajectories of biochemical activity relatively well when their associated points in parameter space are systematically perturbed. This can be inferred by observing the excess of small average Dvalues associated to each reaction parameter. Interestingly, a notable statistical tendency with respect to the system's dynamical behavior was revealed. For example, the signaling network was found to be moderately and extremely sensitive to random perturbations in few reaction parameters. For instance, the parameters related to the Dephosphorylation Rate of the IKKcomplex, and the Maximal Transcriptional rates and Transcriptional Efficiencies associated to the TnfĪ± and Cxcl10 genes can be categorized as moderately sensitive parameters, with average Dvalues ranging between 0.09 and 0.11. On the extreme side of the sensitivity spectrum, we found that the parameters related to the Production and Degradation rates of the TLR4 Susceptible Form, the Association and Dissociation Rates between Phosphorylated IKKcomplex and IkBNFkB, and the Dissociation Rate between IkB and NFkB, represent extremely critical (sensitive) points of the proposed reaction mechanism, with average Dvalues ranging between 0.12 and 0.58. Furthermore, our statistical analysis also revealed that the variability of the Dvalues for those parameters categorized as moderately and extremely sensitive were found to be extremely large, as indicated by both the height of bars and their corresponding whiskers. This result strongly suggests that the robustness properties of the network can be highly variable depending on its current position in biochemical reaction space. It is also interesting to analyze our simulation results from the viewpoint of sloppy and stiff multidimensional parameter spaces [11, 12]. According to this wellsupported theoretical framework, we may conclude that our proposed reaction scheme functions as a highly sloppy information processing system capable of performing robustly, despite undergoing simultaneous random perturbations in its internal reaction parameters. However, some stiff axes were found to be a defining feature of this multidimensional space, along which random perturbations lead predominantly to dramatic changes in the global dynamical behavior of the system. Therefore, such stiff axes in biochemical reaction space constitute key variational constraints of the proposed reaction mechanism. Following this direction, our simulation results strongly suggest that those biochemical processes relying on the reaction parameters identified as critical points of the network, should represent the rate limiting steps that most effectively control the global dynamical behavior of the system. We thus predict that such critical reaction steps represent ideal candidates for manipulating the dynamic activity of the TLR4 signaling network via multitarget therapeutic strategies, which might provide the means for modulating quantitatively innate immune cellular responses in an efficient manner.
Variability of key individual dynamical trajectories
Further statistical analyses were performed to characterize the variability of the dynamical trajectory displayed by each individual reaction species modeled, upon systematic perturbation of the entire biochemical reaction space. We calculated the coefficient of variation of the discrepancies of individual trajectories from the corresponding reference trajectory. A simple Euclidean metric was implemented to evaluate such discrepancies (see the "Models and computational framework" section below); again, this was carried out for each reference trajectory included in the final ensemble of 100 trajectories simulated. In this analysis we focused on those dynamical trajectories categorized as robust/insensitive according to our previous MPSA. This analysis provides primary information on key variational constraints in the network's dynamical behavior. Figure 6 illustrates the results of our statistical analysis, wherein a highly heterogenous spectrum of variability can be readily appreciated, indicating that not all the dynamical trajectories of individual reaction species tend to vary similarly upon global perturbation of the biochemical reaction space. Notably, the temporal trajectory of some reaction species was found to be more variable than others. Note for example that the most upstream (i.e. TLR4 activated form) and the most downstream reaction species (TnfĪ± and Cxcl10) along the signaling cascades modeled, exhibited a remarkable tendecy to vary upon global perturbations. This can be explained by noting that the reaction rules underlying ligandreceptor and transcriptional kinetics involve highly nonlinear processes related to cooperativity interactions and saturation phenomena. Therefore, it should be expected for multiple combination of perturbations in the kinetic parameters driving such non/linear reaction processes to exert drastic changes in the dynamical trajectories of the system. Alternatively, we also found that some intermediate reaction species along the signaling cascades analyzed exhibit a considerable tendency to vary; although some notable differences were observed. For example, a large number of species associated to the MyD88dependent signaling pathway, namely, TABTAK, TABTAKp, P38pn, IKKc, IKKcp, NFĪŗ B, NFĪŗ Bn, TpL2p, were found to vary considerably; whereas only the reaction species TRAM, in the alternative reaction cascade downstream the TLR4, was found to vary significantly. It is tempting to speculate on these results based on the fact that a larger density of coupled biochemical reactions along the MyD88dependent signaling pathway occurs during the time scale considered in our simulations. From this, it then follows that a stronger dynamical coupling of biochemical reactions (functional dependencies/linkages) through this pathway might lead to considerably larger effects when multiple perturbations are propagated dynamically.
Comparison of total parameter variation spectra
Finally, we sought to quantify the capacity of the network of absorbing large fluctuations in internal reaction parameters, and in different regions of biochemical reaction space. We assessed and compared the spectra of total parameter variation (T) for those configurations that were identified as robust and fragile (sensitive) according to our previous MPSA. This measure provides a quantitative notion of the order of magnitude in the variation of a perturbed parameter configuration obtained from a reference one (see the "Models and computational framework" section below). The analysis was performed in each of the 100 dynamical trajectories previously assembled. In Figure 7, every vertical line of points illustrated in each panel stands for a distribution of T values calculated for each reference dynamical trajectory defined by a given point in parameter space. To test for statistical differences between the two spectra shown in Figure 7, we ran MannWhitney tests between robust and fragile distributions. Of the 100 statistical tests performed, we found that 67% of them yielded pvalues < 0.05, thus indicating that, in general, the two spectra tend to differ significantly. However, a simple graphical comparison between the two spectra indicates that a similar global tendency appears to exist (see ranges of variation, for example). In other words, this seems to suggest that the capacity of the signal transduction network of absorbing random perturbations in the whole set of internal reaction parameters may be quite similar under both robust and fragile conditions. At first glance, this observation appears counterintuitive, because it would be expected that for those perturbed configurations categorized as robust/insensitive, small quantitative departures from the reference parameter configuration should be a prevailing statistical regularity. This observation is consistent with the idea that the robust dynamical performance of the network should be more heavily dependent on the direction towards which random perturbations are induced in the biochemical reaction space, than on the magnitude of the perturbation itself.
Robustness of particular inputoutput maps: effects of local and global perturbations at the level of individual transcriptional outputs
Upon extensive exploration and statistical characterization of general robustness properties inferred from hypothetical, but biologically plausible dynamical trajectories displayed by the network, we then focused on a detailed analysis of particular inputoutput maps embedded in the model reaction scheme.
Specifically, we sought to characterize the robustness of the temporal trajectory of the two transcriptional readouts incorporated in the signaling network model. TnfĪ± and Cxcl10 represent crucial outputs required for the appropriate development of proinflammatory responses, which are critically modulated by the upstream reaction cascades activated upon LPS stimulation. Figure 8 shows the temporal profile for the transcriptional activation of these genes. It is worth emphasizing that our modeling framework only accounts for transcriptional activation processes during the short time window simulated. Accordingly, it was assumed that transcriptional activation was mainly driven by the activity of ERK, P38, NFĪŗ B, JNK, and IRF. Consequently, transcriptional repression effects were deliberately neglected. Before conducting the perturbation experiments, several random explorations of the biochemical reaction space were first performed to identify different parameter configurations capable of reproducing the reported experimental data (Figure 8, blackdashed trajectories). Monte Carlo simulations were thus performed taking as reference the set of kinetic data (both initial conditions and parameter ranges) used for simulating the ensemble of reference dynamical trajectories previously constructed (see Additional file 1). Moreover, based on this same reference ensemble of kinetic data, random searches through parameter space via a pseudorandom search algorithm (PRSA) were also performed, as described in [55]. This time, we constructed a small ensemble of 10 reference parameter configurations widely scattered in biochemical reaction space (see Additional file); each reference parameter configuration was able to reproduce relatively well, in quantitative terms, the empirical data (Figure 8, colorcoded trajectories). We focused on this ensemble for conducting local and global perturbation analysis with respect to the transcriptional activation profiles of the TnfĪ± and Cxcl10 proinflammatory genes.
Local perturbation analysis of transcriptional outputs
The computational strategy designed for systematically exploring orthogonal (i.e. local) perturbations in each of the 10 reference points previously sampled is described below (see the "Models and computational framework" section below). In this analysis, perturbations were restricted to those reaction parameters sustaining only signaling fluxes, while transcriptional parameters were maintained unperturbed. In this way, we were able to analyze the quantitative effects at the level of transcriptional readouts of small perturbations impinging upon single reaction kinetics of the upstream signaling cascades. Tables 1 and 2 summarize the calculated overall state sensitivity coefficients for a range of magnitudes of the perturbations induced (ā Ī P ā {0.1, 0.2, 0.3, 0.4, 0.5}, see the "Models and computational framework" section below) in each single reaction parameter. Importantly, the calculated coefficients were found to be remarkably small in comparison to other calculated values reported in a previous simulation study of the MAPK signaling module (see [55]). This observation clearly indicates that the reaction mechanism underlying TLR4mediated signaling behaves as a robust information processing system in the face of small quantitative fluctuations in individual reaction parameters. In other words, this inherent robust condition endows the signaling network, under very particular mutational conditions (i.e point mutations) and within a rather short time window, with the capacity of converting an external stimuli into highly reproducible transcriptional readouts. We note that the two transcriptional readouts considered tend to exhibit the same pattern of sensitivity to most reaction parameters. Interestingly, however, some parameters were found to be more determinant for one transcriptional readout than for the other. For example, the temporal expression profile of Tnfa was found to be relatively more sensitive to variations in k 31f (Association Rate between TLR4I1 and I2), k 35r (Dissociation Rate of the Complex TLR4I1I2I3TT from RIP), ksa (Transition Rate from Susceptible to Active TLR4), and k 14r (Dissociation Rate of the Complex TABTAKp  MKK3/6); whereas the dynamical trajectory of the transcriptional readout of Cxcl10 exhibited a more pronounced sensitivity to perturbations in k 40f (Association Rate of TLR4I1I2I3TTTTBK1 with IR F), k 42f (Dimerization Rate of IRFp), kas (Transition Rate from Active to Susceptible TLR4), and k 1r (Dissociation Rate of MYD88/Mal  TLR4). Although informative, results from local (orthogonal) perturbation analysis provide only limited insight on the variational constraints and systemslevel properties of the signal transduction network. We next performed a global perturbation analysis to this aim, as described below.
Revealing the global perturbation landscapes of transcriptional outputs
In this analysis, we considered global (nonorthogonal) perturbations systematically induced on each of the 10 reference points distributed in biochemical reaction space. We generated 5000 perturbed configurations from each reference point, following our global perturbation strategy. This set of perturbation experiments were carried out to characterize the global perturbation landscapes based on well identified inputoutput relationships. We thus calculated a large number of Dstatistics (obtained via the MPSA approach) for each simulated transcriptional output. That is, Dstatistics were computed for each parameter in 10 different regions of the biochemical reaction space. Briefly, the perturbation landscapes illustrated in Figure 9 provide systemslevel insights into the robust properties and information processing capabilities of the network, but this time in terms of particular inputouput maps embedded in the reaction system. In these landscapes the reaction parameter axis, which ranges between 1108, indicates the set of model parameters influencing directly or indirectly the transcriptional activation of each gene; whereas the parameter configuration axis, ranging between 110, describes the reference ensemble comprising 10 parameter configurations. Remarkable statistical regularities were found when analyzing the perturbation landscapes. 1) In general, the Dvalues associated to many parameters can be highly variable; these were found to be heavily dependent on the parameter configuration tested. This observation indicates that the reproducibility of particular transcriptional readouts in the face of global perturbations is strongly dependent on the region in biochemical reaction space occupied by the signal transduction network. For example, under some parameter regimes the dynamical trajectory of the transcriptional output may be more sensitive to variation in some parameters than in others. 2) In particular, the perturbation landscape for TnfĪ± tend to exhibit and extremely rough topography, with an excess of large Dvalues (> 0.3) distributed heterogenously all over the surface; whereas the landscape associated to Cxcl10 was found to be particularly flat, with just few regions displaying large Dvalues (> 0.3). These findings provide convincing statistical evidence supporting the idea that the transcriptional readout of TnfĪ± should be remarkably more sensitive in the face of global fluctuations in internal reaction parameters than that expected for the transcriptional readout of Cxcl10. Moreover, in the case of the transcriptional readout of TnfĪ±, it is notable the way in which Dvalues associated to a given parameter tend to fluctuate depending on the parameter configuration, that is, the position in parameter space. Alternatively, most Dvalues calculated with respect to the transcriptional readout of Cxcl10 were found to fluctuate only slightly; in general D statistics exhibit a rather invariable tendency across different regions of parameter space. Only in few cases (specific points in biochemical reaction space) considerably large Dvalues were found in the landscape calculated for Cxcl10. For instance, the analysis reveal that only a small fraction of the whole set of reaction parameters appears to most effectively control the transcriptional readout of Cxcl10, which include the following parameters: k 21cat (the Dissociation Rate of IĪŗBNFĪŗB), k 22f (the Import Rate to Nucleus of NFĪŗB), Ī± 2 (Transcriptional Regulatory Strength of IRFpp* over Cxcl10), Ī² 2 (Transcriptional Regulatory Strength of NF Īŗ B over Cxcl10), V A 2 (Cooperativity Effects of IRFpp* on Cxcl10), V B 2 (Cooperativity Effects of NF Īŗ B on Cxcl10), K Ī² 2 (MichaelisMentenConstant Related to Cxcl10 Transcription), k_{ d }Cxcl (Degradation Rate of Cxcl10 mRNA), TmaxCxcl (Max. Transcriptional Rate of Cxcl10), and ĻCxcl (Transcriptional Efficiency of the Cxcl10 Promoter). Taken together, these simulation results point to the idea that the sensitivity/robustness of a given gene expression pattern should be strongly dependent on the architecture of the signaling fluxes influencing directly or indirectly its transcriptional activation. Following this line of arguments, it is interesting to note that the transcriptional activation of TnfĪ±, within our short time scale simulated, relies indirectly on the intranuclear activation of ERK, P38, and JNK, which activates AP1, which in turn activates the transcription of TnfĪ± along with NFĪŗ B; whereas, the transcriptional activation of Cxcl10 only relies on NFĪŗ B and IRF. Under these considerations, it should be clear that the density of signaling fluxes exerting control over the activation of TnfĪ± far exceeds the density of fluxes influencing the activation of Cxcl10. Our observations thus point to the idea that the propagation of multiple perturbations along the reaction cascades should differentially impact the temporal trajectory of the transcriptional readouts of TnfĪ± and Cxcl10. Nevertheless, such apparent differences observed in the topography of the perturbation landscapes are likely to vanish under different molecular scenarios. For example, feedback control or systematically correlated perturbations among subsets of parameters may lead to rather similar perturbation landscapes.
Discussion
The main purpose of this in silico work was to explore whether important systemlevel attributes of a complex biomolecular network were strongly conditioned by the type of signaling tasks (i.e. particular dynamical regimes of molecular activity) simulated. Specifically, our computational approach permitted us an unbiased statistical assessment of the robustness properties, as well as the information processing capabilities, of the canonical reaction mechanism underlying TLR4mediated signal transduction events. This was achieved by considering a broad spectrum of plausible dynamical behaviors displayed by the network (including wild type phenotypes), which are likely encountered in any cell lineage (i.e. macrophage) under diverse physiological conditions. This is the rationale behind our work, and we highlight that these considerations have been largely underappreciated in previous studies of network robustness. Recent investigations, however, have stressed the importance of assessing the spectrum of variational constraints (i.e. robustness, evolvability, epistasis, etc.) of complex developmental regulatory networks under different hypothetical and observable dynamical regimes [13, 56]. Our work thus differs considerably from recent computational studies wherein heavy emphasis have been placed on the characterization of robustness of particular intracellular networks under rather limited biological circumstances [17, 18, 27, 28].
To summarize, our numerical findings strongly suggest that the canonical TLR4 signaling network that drives crucial innate immune cellular responses in macrophages, should be operative in widely scattered regions of the biochemical reaction space; a robust property that allows the network to perform complex signaling tasks in a highly reproducible manner under rather different regimes of molecular activity, and when facing multiple kinetic uncertainties.
Deliberately, we have restricted our model signal transduction network to a simple biochemical reaction mechanism. Importantly, the design principle (topology) of the network was mathematically represented by means of basic reaction schemes defined in terms of mass action law and Hill saturation kinetics. Accordingly, information processing in our model network takes place only through the kinetic coupling of multiple, but rather simple, reaction rules accounting for ligandreceptor interaction, association and dissociation events between single or multiple reaction species, import/export fluxes between cellular compartments, enzymecatalyzed reactions, and transcriptional control. Elaborated regulatory schemes, such as inhibitory reactions or feedback control, were not accounted for in our modeling framework. This is because within our narrow temporal window, in which immediate immune cellular responses are elicited, signal propagation is thought to be controlled in its entirety by the intrinsic crosstalking of MyD88dependent and TRIFdependent reaction cascades (see [46, 47] and references therein). Therefore, within our simulated time window, emphasis was not placed on the complex negative feedback control arising within the NFĪŗ B regulatory module, which is triggered by a wide spectrum of proinflammatory stimuli [57]. The many possible roles of negative feedback control deployed by the NFĪŗ B regulatory module under different cellular contexts have been a central theme of investigation in intracellular signaling ([51, 57]); this issue, however, was beyond the scope of our study. Nevertheless, we acknowledge that our results on the robustness properties and information processing capabilities of the TLR4 signaling network are expected to differ considerably under a different mathematical representation of the reaction topology, wherein positive/negative feedback regulation taking place at any point along the signaling cascade were accounted for. This should come as no surprise, since the crucial role of such elaborated regulatory schemes in any signal transduction system has been well documented (see for example [16, 58], and references therein).
Once having clarified the scope of our study, and more specifically the range of validity of our numerical experiments, we would like to discuss the biological implications of our major findings. Specifically, our global perturbation analyses provide valuable information with respect to plausible variational constraints arising in the system as a result of its canonical design principle. For example, our simulation results indicate the presence of key rate limiting steps that seem to most effectively control the dynamical behavior of the signal transduction network. In particular, statistical analyses from nonorthogonal pertubation experiments clearly show that the global behavior of the system is tightly controlled by only a tiny fraction of the reaction steps embedded in the whole reaction mechanism (see Figure 10). On the other hand, our analyses from global perturbation experiments based on the temporal profiles of transcriptional activation of the two proinflammatory genes modeled (TnfĪ± and Cxcl10) indicate two very distinct scenarios of signaling control. Firstly, the control of the transcriptional readout of TnfĪ± is surprinsingly distributed throughout the whole reaction network (see Figure 9, top panel, and Figure 11). Hence, the transcriptional activation of TnfĪ± should be tightly kinetically involved by virtue of the signaling fluxes displayed by the network upon LPS stimulation. Secondly, the control of the transcriptional readout of Cxcl10 was found to be sparsely distributed in the reaction network, with just few reaction steps critically involved in this signaling ouput (see Figure 9, bottom panel, and Figure 11). Taken together, our results provide mechanistic insights on complex aspects of intracellular signaling in the context of innate immune cellular responses, which might be universal principles of cellular information processing. Overall, our numerical experiments agree well with results from a recently published simulation work [20] indicating that care should be taken when analyzing the robustness properties of any biomolecular network, as these can be heavily dependent on both the quantitative outputs being evaluated, and the current kinetic status of the system (i.e. the position in biochemical reaction space). Finally, and perhaps most importantly, our computational study strongly suggests that the development of effective therapeutic strategies aimed at modulating particular cellular responses, such as metabolic fluxes, signaling and transcriptional outputs, should place heavy emphasis on the architecture of the underlying biomolecular systems. Interestingly, congruent with our findings, accumulating numerical evidence have demostrated that cellular information processing seems to emerge mainly from highly nonlinear dynamics, as well as synergistic/antagonistic interactions among system's components, which can not be resolved by intuitive reasoning alone.
Final remarks
Most systems biology studies centered on the structural and functional organization of highlydimensional biomolecular systems point to the general idea that signal transduction networks should display the inherent capacity of accomplishing specific biological tasks in a robust manner (see [12] and references therein). Robustness seems to be a natural property stemming from the evolved design principle of biomolecular networks [7ā9], which allow them to inhabit sloppy parameter spaces wherein system's behavior turn out to be highly sensitive to variation along a few stiff directions, while being remarkably insensitive to variation along a large number of sloppy axes in parameter space [11, 30]. Notably, accurate computational reconstructions of experimentally reported dynamical behaviors of many signal transduction networks have been successfully achieved [20, 51, 55, 57]. Interestingly, standard mathematical representations of the reaction topology of most signaling network models are typically founded on highly nonlinear, but relatively simple, biochemical reaction rules, which despite being an abvious simplification of the underlying biochemistry have proven successful at providing mechanistic insight [20, 51, 55, 57]. This is an intriguing observation from an evolutionary standpoint. This suggests, for example, that the underlying mathematical structure of most signal transduction networks that has been favored over evolution to process in an efficient and robust manner the biochemical information arising in the cell, might simply rely on basic dynamic rules ([59, 60]). Intuitively, the most variable component affecting the temporal variation in the activity of the molecular species involved in a certain signaling event would be the number of the contributing reaction velocities to a particular flux. Following this line of arguments, it is tempting to speculate on the possibility that the deterministic component of the dynamical trajectories displayed by most signaling networks might have been the result of selection for simple biochemical reaction rules built, for example, upon mass action and Hilllike saturation kinetics.
Models and computational framework
The mathematical representation of the canonical reaction network retrieved from the literature, and the whole set of numerical experiments that are described below were implemented in Mathematica^{Ā®} 6.0.
Mathematical formulation of the signal transduction network in the language of dynamical systems
A signal transduction network can be appropriately conceived in dynamical terms, whose internal regulatory schemes, reaction rules and associated control parameters underlying the trajectories of the system can be formulated, as a first approximation, via basic principles of biochemical reaction. Ordinary Differential Equationsbased models grounded on mass action law (first and second reaction kinetic orders) and Hill saturation kinetics, provide a suitable macroscopic approximation to intracellular signal transduction dynamics (fluxes), as well as transcriptional phenomena. Under this modeling formalism, a biochemical reaction network can be described on the basis of the state space formulation with the following mathematical constructs:

1
A network with n reaction species is represented by the state vector:
(1)In our case, the TLR4 signaling network model incorporates n = 76 reaction species, including receptors, adapters, kinases, transcription factors and mRNAs.

2
The biochemical reaction parameters controlling the signaling flux through the network and transcriptional processes of target genes are incorporated in the reaction vector
(2)With m = 116 internal reaction coefficients for the TLR4 network. Here, Ī encompasses a wide spectrum of parameters of different biochemical nature, ranging from transition rates between receptor states (susceptible ā activated), production and degradation rates of receptors, association/dissociation rates among intracelular molecular species, phospho/dephosphorylation rates, nuclear import/export rates, maximal transcriptional rates, transcriptional efficiencies, MichaelesMenten constants, cooperative coefficients, and mRNA degradation rates (see Additional file 1 for a detailed description of these parameters and their assigned range of values).
According to the above mathematical expressions, the state space nonlinear representation of the biochemical reaction network is thus given by:
Where f(Ā·) defines a nonlinear state transition function accounting for reaction velocities or fluxes (see Additional file 1 for a detailed description of the dynamical system), which can be grounded on mass action laws and/or Hill kinetics, according to the reaction mechanisms modeled (receptor activation kinetics, binding and enzymatic reactions, or transcriptional dynamics). Y_{0} represents the vector of initial concentrations for the reaction species at time t_{0}; g(Ā·) defines a measurement function which is solved numerically; whereas X ā ā^{n}gives the measurement output vector representing the concentration of the reaction species at a given point in time. In our case, the TLR4 signaling network model accounts for 76 dynamic variables (reaction species, Y_{ j }ā j ā (1, 2, 3, ....,76)), 32 of which were assigned zero initial conditions (see Additional file 1).
Multiparametric sensitivity analysis (MPSA): a combination of uncertainty and sensitivity analyses
The MPSA approach was first introduced by Hornberger and Chang [61, 62]. This is a computational strategy specially suitable for characterizing the relative importance of the parameters of a multidimensional mathematical model. Additionally, a MPSA also provides the means for the identification of possible correlations within the system under study. The basic idea of the MPSA method is to generate ensembles of different input parameter configurations for systematically evaluating the range of responses (outputs) accessible to a given mathematical model. This can be achieved via Monte Carlo simulations, in which the model is run iteratively using sets of parameters drawn randomly from prespecified distributions. Since the natural distributions of parameter values for actual biological networks are completely unknown, we thus implemented uniform probability distributions in our simulations. Since our biochemical reaction network inhabits a multidimensional parameter space encompassing 116 biochemical reaction axes, which needs to be sampled efficiently, we thus implemented a Latin Hypercube Sampling (LHS) method so as to generate representative ensembles of different parameter configurations in both the immediate and distant vicinity of a reference point in parameter space. This initial approach is aimed at injecting uncertainty into the model's inputs, and has thus been coined the term Uncertainty Analysis [23]. We then made use of these ensembles of parameter configurations obtained from previously assembled reference points in parameter space, to statistically characterize the robust properties of the canonical reaction topology underlying TLR4mediated signaling processes in macrophages; a computational approach known as Sensitivity Analysis [23ā26]. We define a decision rule that allowed for the classification of each model evaluation as either acceptable (robust) or unacceptable (nonrobust). Further, statistical analyses, based on KolmogorovSmirnov (KS) tests, were performed on the occurences of acceptable and unacceptable cases, which are summarized for each model parameter. The larger the difference between the cumulative frequencies of the two cases, which is reflected by large values of the D statistic obtained from KS tests, the more significant is a given parameter. The MPSA method can be summarized in the following steps:

1
Selection of reference parameter configurations (vectors) to be perturbed.

2
Set a relatively large range of variation for each model parameter in order to account for a wide spectrum of biologically plausible perturbations (i.e. single or combined mutations, thermal fluctuations, etc.). In our case, a perturbation variable, Ļ, was sampled in this way: Ļ ~ U(1, 1); a perturbation function was then applied over a reference parameter i, Īø_{i, ref}, in order to obtain a newly perturbed parameter Īø_{i, pert}= 10^{Ļ}* Īø_{i, ref}

3
Under this perturbation strategy we initially generated seeds of LHS matrices with the basic structure as illustrated in Table 3, wherein a row vector stands for a perturbed parameter configuration. In this way, the matrix shown represents a set of input vectors distributed in parameter space in the vicinity (either immediate or distant) of a previously defined reference point; this matrix is assembled via the LHS strategy, and is designed for the systematic evaluation of the model's output. Note that in our case 5000 parameter configurations were generated from a previously defined reference point in parameter space. We constructed ensembles of 100 LHS matrices for evaluating the robust information processing capabilities of the network in different regions of biochemical reaction space. We also constructed an ensemble of 10 LHS matrices in order to test for the robustness properties of the two experimentally reported transcriptional outputs of the network. Before simulating each previously assembled seed LHS matrix, however, we permuted the elements of each column of a matrix as illustrated in Table 4. In this way, permuting the matrices permitted us to avoid any kind of bias in model evaluations.

4
Each LHS matrix was then simulated, and the corresponding discrepancy function evaluated, which is of the form:
(5)With J denoting a reference dynamical trajectory with associated parameter configuration Ī_{ ref }, and Ī_{ pert }being a perturbed version of it obtained from a LHS matrix; with h ā (1, ..., 5000) being any evaluation of the discrepancy function.

5
This step implied determining whether a given perturbed parameter configuration was acceptable (robust) or unacceptable (sensitive) by comparing the discrepancy function value to a given threshold. If the discrepancy function value was found to be below the threshold value the evaluated parameter configuration was then classified as acceptable; otherwhise it was classified as unacceptable. Previous computational works suggest that results from MPSA should not be affected considerably with the choice of a given discrepancy function [25, 53, 63]. Here we implemented the average of the discrepancy function as our threshold value, defined as follows:
(6) 
6
After systematic evaluation of the discrepancy function for each LHS matrix, statistical assessment followed in order to determine whether a given parameter configuration was deemed either acceptable or unacceptable. To do this, we applied the KS test to assess the global sensitivity of the system's output with respect to perturbations targeting individual parameters. The KS test provides the means for evaluating the cumulative frequency of the observations (parameter values) as a function of class, and calculate the maximum vertical distance between cumulative frequency distribution curves for m acceptable and n unacceptable cases of any given parameter Īø_{ j }. This is obtained by calculating the D statistic, which is defined in this way:
(7)
Where S(Īø_{ j }) and represent the cumulative frequency functions corresponding to acceptable and unacceptable cases, respectively, with Īø_{ j }being any reaction parameter of the signal transduction network. Importantly, this estimator provides a robust quantitative notion for the sensitivity/robustness of the network model to random perturbations of the reaction parameters. The higher the Dvalue the more sensitive is the dynamical behavior of the network model with respect to variation of a given parameter, when the remaining parameters (i.e. biochemical background of the network) are also varied.
Total parameter variation
The total parameter variation estimator provides a quantitative notion of the order of magnitude in the variation of a perturbed parameter configuration obtained from a reference one. This estimator is defined as:
We calculated all T values for those perturbed parameter configurations that were deemed either robust or fragile obtained from a reference parameter configuration, Ī_{ ref }, exhibiting a global dynamical trajectory J (ā J ā (1, ......., 100))
Local and global perturbation analysis of inputoutput maps
Here, it is described the global (nonorthogonal) and local (orthogonal) perturbation strategies with respect to particular inputs of the network. Specifically, we assessed the reproducibility of the network's transcriptional outputs (TnfĪ± and Cxcl10) in the face of single (local) and multiple simulataneous (global) perturbations targeting the internal reaction parameters of the network. Local perturbations were quantified via the following metric:
This can be referred to as an "Overal Senstitivity Index" [17, 55], being s_{ i }(t) defined as follows:
With the variable Output representing either TnfĪ± or Cxcl10, and Īø_{ j }giving any internal reaction parameter j belonging to a reference parameter configuration Ī_{ ref }. In this local perturbation experiment we focused only on those reaction parameters accounting for signaling fluxes from the cell membrane surface down to the nucleus; that is to say, perturbations targeting transcriptional parameters were not simulated. On the other hand, to perform global perturbation analysis based on the effects of multiple induced perturbations at the level of transcriptional readouts, we implemented the same MPSA described above. However, in this case we followed a different perturbation scheme. We chose a perturbation variable, Ļ, to be sampled in this way: Ļ ~ U( 1, 1); a perturbation function was then applied over a reference parameter i, Īø_{i, ref}, in order to obtain a newly perturbed parameter Īø_{i, pert}= 5^{Ļ}* Īø_{i, ref}. In this experiment, we perturbed the entire biochemical reaction space, including internal reaction parameters related to both signaling fluxes and transcriptional processes, and assessed the effects of the perturbations at the level of each transcriptional readout separately (TnfĪ± or Cxcl10), via the following discrepancy function:
This time, our criterion used to categorize a parameter configuration as acceptable or unacceptable was based on the following threshold value:
Given in arbitrary units of discrepancy, this threshold was selected upon a detailed analysis of both the qualitative and quantitative effects of the perturbations on the temporal dynamics of the transcriptional readouts.
References
Bhalla US, Iyengar R: Emergent Properties of Networks of Biological Signaling Pathways. Science. 1999, 283: 381387. 10.1126/science.283.5400.381.
Weng G, Bhalla US, Iyengar R: Complexity in Biological Signaling Systems. Science. 1999, 284: 9296. 10.1126/science.284.5411.92.
Kholodenko B: Cell Signaling Dynamics in Time and Space. Nat Rev Mol Cell Biol. 2006, 281 (29): 1992519938.
Hlavacek WS, Faeder JR: The Complexity of Cell Signaling and the Need for a New Mechanics. Science Signaling. 2009, 2 (81): pe4610.1126/scisignal.281pe46.
Klemm K, Bornholdt S: Topology of Biological Networks and Reliability of Information Processing. PNAS. 2005, 102 (51): 1841418419. 10.1073/pnas.0509132102.
Helikar T, Konvalina J, Heidel J, Rogers JA: Emergent DecisionMaking in Biological Signal Transduction Networks. PNAS. 2008, 105 (6): 19131918. 10.1073/pnas.0705088105.
Alon U, Surette MG, Barkai N, Leibler S: Robustness in Bacterial Chemotaxis. Nature. 1999, 397 (6715): 168171. 10.1038/16483.
von Dassow G, Meir E, Munro EM, Odell GM: The Segment Polarity Network is a Robust Developmental Module. Nature. 2000, 406: 188192. 10.1038/35018085.
Meir E, von Dassow G, Munro E, Odel GM: Robustness, Flexibility, and the Role of Lateral Inhibition in the Neurogenic Network. Curr Biol. 2002, 12: 778786. 10.1016/S09609822(02)008394.
Pribyl M, Muratov CB, Shvartsman SY: Transitions in the Model of Epithelial Patterning. Dev Dyn. 2003, 226 (1): 155159. 10.1002/dvdy.10218.
Daniels BC, Chen YJ, Sethna JP, Gutenkunst RN, Myers CR: Sloppiness, Robustness, and Evolvability in Systems Biology. Curr Opin Biotech. 2008, 19: 1:7
Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP: Universally Sloppy Parameter Sensitivities in Systems Biology Models. PLoS Comput Biol. 2007, 3 (10): e18910.1371/journal.pcbi.0030189.
GutiĆ©rrez J: A Developmental Systems Perspective on Epistasis: Computational Exploration of Mutational Interactions in Model Developmental Regulatory Networks. PLoS ONE. 2009, 4 (9): e682310.1371/journal.pone.0006823.
Hornberg JJ, Binder B, Bruggeman FJ, Schoeberl B: Control of MAPK Signaling: From Complexity to what Really Matters. Oncogene. 2005, 24: 55335542. 10.1038/sj.onc.1208817.
Heinrich R, Neel BG, Rapoport TA: Mathematical Models of Protein Kinase Signal Transduction. Molecular Cell. 2002, 9: 957970. 10.1016/S10972765(02)005282.
Stelling J, Sauer U, Szallasi Z, Doyle FJI, Doyle J: Robustness of Cellular Functions. Cell. 2004, 118 (6): 675685. 10.1016/j.cell.2004.09.008.
Stelling J, Gilles ED, Doyle FJ: Robustness Properties of Circadian Clock Architectures. PNAS. 2004, 101 (36): 1321013215. 10.1073/pnas.0401463101.
Kurata H, Tanaka T, Ohnishi F: Mathematical Identification of Critical Reactions in the Interlocked Feedback Model. PLoS ONE. 2007, 2 (10): e110310.1371/journal.pone.0001103.
Zou X, Liu M, Pan Z: Robustness Analysis of EGFR Signaling Network with a MultiObjective Evolutionary Algorithm. BioSystems. 2008, 91: 245261. 10.1016/j.biosystems.2007.10.001.
Chen WW, Schoeberl B, Jasper PJ, Niepel M: InputOutput Behavior of ErbB Signaling Pathways as Revealed by a Mass Action Model Trained against Dynamic Data. Mol Syst Biol. 2009, 5: 239
Chaves M, Sengupta A, Sontag ED: Geometry and Topology of Parameter Space: Investigating Measures of Robustness in Regulatory Networks. J Math Biol. 2008
Dayarian A, Chaves M, Sontag ED, Sengupta AM: Shape, Size and Robustness: Feasible Regions in the Parameter Space of Biochemical Networks. PLoS Comput Biol. 2009, 5 (1): e100025610.1371/journal.pcbi.1000256.
Marino S, Hogue IB, Ray CJ, Kirschner DE: A Methodology for Performing Global Uncertainty and Sensitivity Analysis in Systems Biology. J Theor Biol. 2008, 254: 178196. 10.1016/j.jtbi.2008.04.011.
Yue H, Brown M, Knowles J, Wang H: Insights into the Behavior of Systems Biology Models from dynamic Sensitivity and Identifiability Analysis: A Case Study of an NFkB Signaling. Molecular BioSystems. 2006, 640: 640649. 10.1039/b609442b.
Cho KH, Shin SY, Kolch W, Wolkenhauer O: Experimental Design in Systems Biology, Based on Parameter Sensitivity Analysis Using a Monte Carlo Method: A Case Study for the TNFalphaMediated NFkB Signal Transduction Pathway. Simulation. 2003, 79 (1112): 115.
Van Riel NAW: Dynamic Modelling and Analysis of Biochemical Networks: Mechanismbased Models and Modelbased Experiments. Briefings in Bioinformatics. 2006, 7 (4): 364374. 10.1093/bib/bbl040.
Hafner M, Koeppl H, Hasler M, Wagner A: Glocal Robustness Analysis and Model Discrimination for Circadian Oscillators. PLoS Comput Biol. 2009, 5 (10): e100053410.1371/journal.pcbi.1000534.
Yoon J, Deisboeck TS: Investigating Differential Dynamics of the MAPK Signaling Cascade Using a MultiParametric Global Sensitivity analysis. PLoS ONE. 2009, 4 (2): e456010.1371/journal.pone.0004560.
Nijhout HF, Berg AM, Gibson WT: A Mechanistic Study of Evolvability Using the MitogenActivated Protein Kinase Cascade. Evol Devel. 2003, 5 (3): 281294. 10.1046/j.1525142X.2003.03035.x.
Brown KS, Hill CC, Calero GA, Myers CR, Lee KH: The Statistical Mechanics of Complex Signaling Networks: Nerve Growth Factor Signaling. Phys Biol. 2004, 1: 184195. 10.1088/14783967/1/3/006.
Banerjee A, Gerondakis S: Coordinating TLRActivated Signaling Pathways in Cells of the Immune System. Immunol Cell Biol. 2007, 85: 420424. 10.1038/sj.icb.7100098.
Lu YC, Yeh WC, Ohashi PS: LPS/TLR4 signal Transduction Network. Cytokine. 2008, 42: 145151. 10.1016/j.cyto.2008.01.006.
Krishnan J, Selvarajoo K, Tsuchiya M, Lee G, Choi S: Tolllike Receptor Signal Transduction. Experimental and Molecular Medicine. 2007, 39 (4): 421438.
Oda K, Kitano H: A Comprehensive Map of the Tolllike Receptor Signaling Network. Mol Syst Biol. 2006,2006.0015,
Kaway T, Akira S: The Roles of TLRs, RLRs and NLRs in Pathogen Recognition. Int Immunol. 2009, 21 (4): 317337. 10.1093/intimm/dxp017.
Trinchieri G, Sher A: Cooperation of Tolllike Receptor Signals in Innate Immune Defence. Nat Rev Immunol. 2007, 7 (3): 179190. 10.1038/nri2038.
Shizuo A, Takeda K: Tolllike Receptor Signaling. Nat Rev Immunol. 2004, 4 (7): 499511. 10.1038/nri1391.
Behar M, Dohlman HG, Elston TC: Kinetic Insulation as an Effective Mechanism for Achieving Pathway Specificity in Intracellular Signaling Networks. PNAS. 2007, 104 (41): 1614616151. 10.1073/pnas.0703894104.
Shankaran H, Wiley HS, Resat H: Receptor Downregulation and Desensitization Enhance the Information Processing Ability of Signaling Receptors. BMC Syst Biol. 2007, 1: 4810.1186/17520509148.
Shankaran H, Resat H, Wiley HS: Cell Surface Receptors for Signal Transduction and Ligand Transport: A Design Principles Study. PLoS Comput Biol. 2007, 3 (6): e10110.1371/journal.pcbi.0030101.
Gunawardena J: Signals and Systems: Towards a Systems Biology of Signal Transduction. Proceedings of the IEEE. 2008, 96 (8): 13861397. 10.1109/JPROC.2008.925413.
Natarajan M, Lin KM, Hsueh RC, Sternweis PC, Ranganathan R: A Global Analysis of CrossTalk in a Mammalian Cellular Signalling Network. Nat Cell Biol. 2006, 8 (6): 571580. 10.1038/ncb1418.
Hsueh RC, Natarajan M, Fraser I, Pond B, Liu J, Mumby S: Deciphering Signaling Outcomes from a System of Complex Networks. Science Signaling. 2009, 2 (71): ra2210.1126/scisignal.2000054.
LĆ¼dtke N, Panzeri S, Brown M, Broomhead DS, Knowles J, Montemurro MA, Kell DB: InformatioTheoretic Sensitivity Analysis: A General Method for Credit Assigment in Complex Networks. J R Soc Interface. 2008, 5: 223235. 10.1098/rsif.2007.1079.
Akira S, Takeda K: Tolllike Receptor Signaling. Nat Rev Immunol. 2004, 4: 499511. 10.1038/nri1391.
Selvarajoo K: Discovering Differential Activation Machinery of the Tolllike Receptor 4 Signaling Pathways in MYD88 Knockouts. FEBS Letters. 2006, 580: 14571464. 10.1016/j.febslet.2006.01.046.
Selvarajoo K, Takada Y, Gohda J, Helmy M: Signaling Flux Redistribution at Tolllike Receptor Pathway Junctions. PLoS ONE. 2008, 3 (10): e343010.1371/journal.pone.0003430.
Riviere B, Ephsteyn Y, Swigon D, Vodovotz Y: A Simple Mathematical Model of Signaling Resulting from the Binding of Lipopolysaccharide with Tolllike Receptor 4 Demonstrates Inherent Precoditioning Behavior. Mathematical Biosciences. 2009, 217: 1926. 10.1016/j.mbs.2008.10.002.
An GC, Faeder JR: Detailed Qualitative Dynamic Knowledge Representation Using BioNetGen Model of TLR4 Signaling and Preconditioning. Mathematical Biosciences. 2009, 217: 5363. 10.1016/j.mbs.2008.08.013.
Renner F, Schmitz ML: Autoregulatory Feedback Loops Terminating the NFĪŗB Response. Trends in Biochemical Sciences. 2009, 34 (3): 128135. 10.1016/j.tibs.2008.12.003.
Covert MW, Leung TH, Gaston JE, Baltimore D: Achieving Stability of LipopolysaccharideInduced NFĪŗB Activation. Science. 2005, 309: 18541857. 10.1126/science.1112304.
Stout RD, Suttles J: Functional Plasticity of Macrophages: Reversible Adaptation to Changing Microenvironments. J Leukoc Biol. 2004, 76 (3): 509513. 10.1189/jlb.0504272.
Zi Z, Cho HC, Sung MH, Xia X: In silico Identification of the key Components and Steps in IFNĪ³ Induced JAKSTAT Signaling Pathway. FEBS Letters. 2005, 579: 11011108. 10.1016/j.febslet.2005.01.009.
Jia J, Yue H, Liu T, Wang H: Global Sensitivity Analysis of Cell Signaling Transduction Networks Based on Latin Hypercube Sampling Method. Bioinformatics and Biomedical Engineering, IEEE Xplore. 2007, 434437. full_text.
Shin SY, Rath O, Choo SM, Fee F: Positive and NegativeFeedback Regulations Coordinate the Dynamic Behavior of the RasRafMEKERK Signal Transduction Pathway. J Cell Sci. 2009, 122: 425435. 10.1242/jcs.036319.
Giurumescu CA, Sternberg PW, Asthagiri AR: Predicting Phenotypic Diversity and the Underlying Quantitative Molecular Transitions. PLoS Comput Biol. 2009, 5 (4): e100035410.1371/journal.pcbi.1000354.
Shannon LW, Barken D, Hoffmann A: Stimulus Specificity of Gene Expression Programs Determined by Temporal Control of IKK Activity. Science. 2005, 309: 18571861. 10.1126/science.1113319.
Kholodenko BN: CellSignaling Dynamics in Time and Space. Nat Rev Mol Cell Biology. 2006, 7: 165176. 10.1038/nrm1838.
Bornholdt S: Systems Biology. Less is More in Modeling Large Genetic Networks. Science. 2005, 310 (5747): 449451. 10.1126/science.1119959.
Selvarajoo K, Tomita M, Tsuchiya M: Can Complex Cellular Processes be Governed by Simple Linear Rules?. J Bioinf Comput Biol. 2009, 7 (1): 243268. 10.1142/S0219720009003947.
Hornberger G, Spear R: An Approach to the Preliminary Analysis of Environmental Systems. J Environ Manage. 2005, 12: 718.
Chang FJ, Delleur JW: Systematic Parameter Estimation of Wathershed Acidification Model. Hydrol Processes. 1992, 6: 2944. 10.1002/hyp.3360060104.
Choi J, Hulseapple SM, Conklin MH, Harvery JW: Modeling CO_{2} Degassing and pH in a StreamAquifer System. J Hydrol. 1998, 209: 297310. 10.1016/S00221694(98)000936.
Acknowledgements
We are greatly indebted to Drs. Masa Tsuchiya and Kumar Selvarajoo for helpful disccusions and critical comments on earlier versions of the manuscript. We also thank anonymous reviewers for insightful comments and helpful suggestions. Special thanks to Koichi Matsuo for kindly providing experimental data on gene expression. JG wishes to acknowledge Grupo dĆ© FĆsica y AstrofĆsica Computacional (FACom) for computing facilities, and Boris A. RodrĆguez for important suggestions and constant technical support. This work was funded in part by Grupo de InmunovirologĆa, SIU, Universidad de Antioquia.
Author information
Authors and Affiliations
Corresponding author
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
JG designed the study, developed the computational framework, performed simulations and statistical analysis, discussed the results and drafted the manuscript. GSL and SUI discussed the results and helped to draft the manuscript. All authors read and approved the final version of the manuscript.
Electronic supplementary material
12976_2009_220_MOESM1_ESM.PDF
Additional file 1: Mathematical structure of the signal transduction network: kinetic parameters, initial conditions, and rate equations. This file contains a detailed description of our modeling framework. Ranges of values for kinetic parameters and initial conditions are given, which were selected according to several computational strategies described in the main text. Our system of rate equations implemented for simulating intracellular fluxes and propagation of kinetic uncertainties through the TLR4 signal transduction network, is also described. (PDF 122 KB)
Authorsā original submitted files for images
Below are the links to the authorsā original submitted files for images.
Rights and permissions
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
About this article
Cite this article
GutiĆ©rrez, J., St Laurent, G. & UrcuquiInchima, S. Propagation of kinetic uncertainties through a canonical topology of the TLR4 signaling network in different regions of biochemical reaction space. Theor Biol Med Model 7, 7 (2010). https://doi.org/10.1186/1742468277
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1742468277