Skip to main content
Advertisement
  • Loading metrics

Multiplexing information flow through dynamic signalling systems

  • Giorgos Minas,

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing

    Affiliation School of Mathematics and Statistics, University of St Andrews, St Andrews, Scotland, United Kingdom

  • Dan J. Woodcock,

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Writing – review & editing

    Affiliation Big Data Institute, University of Oxford, Oxford, England, United Kingdom

  • Louise Ashall,

    Roles Conceptualization, Data curation, Investigation

    Affiliation Systems Microscopy Centre, Division of Molecular and Cellular Function, School of Biology, Faculty of Biology, Medicine and Health, Manchester Academic Health Sciences Centre, University of Manchester, Manchester, England, United Kingdom

  • Claire V. Harper,

    Roles Conceptualization, Data curation, Investigation, Writing – review & editing

    Affiliation Systems Microscopy Centre, Division of Molecular and Cellular Function, School of Biology, Faculty of Biology, Medicine and Health, Manchester Academic Health Sciences Centre, University of Manchester, Manchester, England, United Kingdom

  • Michael R. H. White,

    Roles Conceptualization, Funding acquisition, Investigation, Methodology, Supervision, Writing – review & editing

    Affiliation Systems Microscopy Centre, Division of Molecular and Cellular Function, School of Biology, Faculty of Biology, Medicine and Health, Manchester Academic Health Sciences Centre, University of Manchester, Manchester, England, United Kingdom

  • David A. Rand

    Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Supervision, Validation, Writing – original draft, Writing – review & editing

    D.A.Rand@warwick.ac.uk

    Affiliation Mathematics Institute, University of Warwick, Coventry, England, United Kingdom

Abstract

We consider how a signalling system can act as an information hub by multiplexing information arising from multiple signals. We formally define multiplexing, mathematically characterise which systems can multiplex and how well they can do it. While the results of this paper are theoretical, to motivate the idea of multiplexing, we provide experimental evidence that tentatively suggests that the NF-κB transcription factor can multiplex information about changes in multiple signals. We believe that our theoretical results may resolve the apparent paradox of how a system like NF-κB that regulates cell fate and inflammatory signalling in response to diverse stimuli can appear to have the low information carrying capacity suggested by recent studies on scalar signals. In carrying out our study, we introduce new methods for the analysis of large, nonlinear stochastic dynamic models, and develop computational algorithms that facilitate the calculation of fundamental constructs of information theory such as Kullback–Leibler divergences and sensitivity matrices, and link these methods to a new theory about multiplexing information. We show that many current models such as those of the NF-κB system cannot multiplex effectively and provide models that overcome this limitation using post-transcriptional modifications.

Author summary

Cells use signalling systems to pass on information arising from their ever-changing environment to their processing units. These biochemical networks regulate the transmission of multiple signals within the noisy and complex cellular environment, controlling whether to turn on or off processes of cell defence, death, division, and others. The question of how they actually achieve that becomes particularly critical given that many diseases occur when signalling systems malfunction. In this paper, we develop methodology and computational tools for simulating, measuring and analysing the ability of signalling systems to transmit multi-dimensional signals. We specifically focus on the capacity of signalling systems to simultaneously transmit multiple signals, such as temperature changes, presence and concentration of cytokines, viral and bacterial pathogens or drugs, through a single noisy, dynamic signalling system. We argue that a signalling system can act as an information hub, sending information in a multiplexed fashion rather similar to the way in which telecommunications networks send multiple signals over a shared medium by combining them into one.

Introduction

Signalling systems provide a very important example of cellular information systems since they transmit information arising from inside and outside the cell to the cell’s processing units. For example, it is generally believed that the nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) system uses the information from a large number of input signals (see Fig 1(a)) to regulate gene transcription of more than 500 genes in a highly versatile way [1, 2]. NF-κB regulates cell fate and inflammatory signalling in response to diverse stimuli, including changes in temperature [3], viral and bacterial pathogens, free radicals, cytokines, and growth factors [2]. Thus, we have a situation where both the input signal that encodes information about the cell’s environment, and the gene response are multi-dimensional. Such a system is often referred to as an information hub.

thumbnail
Fig 1. Multiplexing signals through signalling systems.

(a) Cells constantly receive a multitude of different signals in which signalling systems respond by (directly or indirectly) modulating the expression of a number of target genes. These target genes activate or not various pathways of the cell which leads to completely different cell outcomes from cell survival to apoptosis or mitosis. In order for this decision making to be reliable and robust, signalling systems need to have the capacity to multiplex a variety of simultaneously arising signals. (b) Multiplexing is defined as the ability of the signalling system response to identify which of the input signals have changed. In broad terms, strong multiplexing is evident by the probability distributions of the signalling system response in a population of single cells being significantly different for the different regimes of the multi-dimensional signal. On the contrary, poor multiplexing leads to response distribution that are very similar for different signals.

https://doi.org/10.1371/journal.pcbi.1008076.g001

This raises the question of to what extent this process is mediated through the signalling system itself which may have a single transcription factor, rather than through multiple other parallel pathways that also provide information to the genome. Can such a signalling system on its own effectively regulate a relationship between multidimensional inputs and responses that can robustly and reliably modulate decision-making of the claimed versatility without using other pathways? This is the central question we consider here.

We consider this question in terms of multiplexing which we define as follows. We suppose that our system has multiple input signals S1, …Ss and consider how the system responds to changes in them. These signals might, for instance, be changes in temperature or other physical parameters (e.g. pressure or humidity), changes in the level and/or timing pattern of an activator (e.g. tumor necrosis factor-α (TNFα), interleukin 1β (IL-1β), and Lipopolysaccharides (LPS) for NF-κB), and/or drug treatments (e.g. Diclofenac for NF-κB). We say that a signalling system can multiplex these input signals S1, …Ss if one can reliably determine which of these input signals have changed using only the multidimensional response of the target genes (see Fig 1(B)). It is this response that will regulate downstream responses of the cell and therefore the multiplexing capacity is directly measuring a key aspect of how effectively the cell can respond to the multiple inputs.

So far as we are aware this is a new approach and consequently an immediate question is whether there is any experimental evidence that this is the case. We address this below and present some tentative evidence for it. This is useful because it provides useful context to our discussion and gives some helpful insights but we must emphasise that this evidence is very far from proving such a point even though it is highly suggestive.

To address the question of what aspects of the system enable such multiplexing we will introduce a quantity, called the multiplexing capacity, which measures the ability of a noisy signalling system to multiplex a set of signals. Using this we demonstrate that while current models cannot multiplex effectively, biologically natural modifications of them can. We believe this indicates general principles behind biological design.

The underlying concepts that we use are connected to important tools involved in systems identification, sensitivity analysis and information geometry such as the Fisher Information Matrix (FIM) and the Kullback-Leibler divergence [4] and we introduce and calculate a new, but related, sensitivity matrix s that characterises the multiplexing capacity. While sensitivity analysis [57] is an extensive area for deterministic dynamics (e.g. [8, 9]) it is less well developed for stochastic systems (e.g. [1016]). Since it is crucial in our discussion that we take account of realistic levels of stochasticity for the systems we consider, calculating the relevant quantities for complex high-dimensional stochastic systems such as those considered here is therefore a significant mathematical challenge. To overcome this, in the numerical computations we use the pcLNA method [17] that allows fast and accurate computation of key information theoretic quantities, such as Kullback-Leibler divergences and the Fisher Information matrix, for stochastic dynamical systems.

Given that NF-κB has complex oscillatory dynamics, an obvious hypothesis is that it is this dynamical behavior of the system that allows it to act as an information hub. However, we will use our theoretical tools to provide evidence that this is not the case and show that the NF-κB system described by current models cannot multiplex effectively even though it has complex oscillatory dynamics. On the other hand, we will demonstrate how to modify a stochastic model of NF-κB so as to overcome this inability to multiplex. In particular, we show that additional regulated states of NF-κB, which might include differential post-translational modifications and/or differential hetero- and homo-dimerisation, can enable such multiplexing and that the oscillatory dynamics can greatly enrich the multiplexing capacity in this modified model.

Recent important papers studied the information flow through biochemical systems such as the nuclear factor-κB (NF-κB), calcium (Ca2+), and extracellular signal-regulated kinase (ERK). The focus was on measuring how much information is being carried by the signalling systems in terms of the mutual information I(S, R) and the capacity of the channel SP(R|S) [4] where S is the input signal and P(R|S) the probability distribution of a response R. For example, for NF-κB the signal S was the level of TNFα stimulation and R was the level of transcription factor in the nucleus at one or more timepoints. In summary, the channel capacity was estimated to be around 1 bit for static scalar observations in response to one-dimensional stimuli [1822], about 1.5 bits when the dynamical behaviour of the system response is considered [21] and up to 1.7 bits when cell-to-cell heterogeneity is accounted for [23]. A number of recent studies support this core observation and report similar low channel capacities [2429]. The stochastic models that we use reproduce these relatively low levels of mutual information between TNFα level and total transcription factor abundance and also agree with that of the stochastic model in [30]. On the other hand, the modified versions allow significantly greater mutual information for multidimensional inputs.

The logic of our discussion is as follows: Firstly, we discuss some tentative experimental evidence suggesting that NF-κB does multiplex. Then we introduce a method for quantifying the ability of a given stochastic model to multiplex and show that current models are poor at this. We suggest how one can modify the models so as to enable better multiplexing and relate this to known mechanisms in signalling systems. In particular, we provide a modified model that is able to reproduce the behaviour of the simpler of the two multiplexing experimental systems we discuss. We finally discuss how these results relate to the low information capacity found in previous studies. Our analysis gives important insight into how multiplexing can work in a signalling system, however we are not claiming that our model is a true representation of the real biology of those systems.

Does NF-κB multiplex?

To illustrate the above characterisation of multiplexing we consider some experimental evidence. We ask if, by monitoring the response of a set of genes that are direct NF-κB targets (Sect. 5 in S1 Appendix), we can reliably determine the state of a multidimensional input signal. We firstly consider the response of three important genes, EGR1, COX-2 (PTGS2) and IL-8 (CXCL-8), to pulses of varying length, repeated every 100 minutes at two temperatures, 37°C and 40°C and ask if, from the response of these genes, we can determine the temperature and pulsing length. EGR1 regulates the response to growth factors, DNA damage, and ischemia, preventing tumor formation by activating p53/TP53 and TGFB1. COX-2 is responsible for production of inflammatory prosta-glandins. We include the chemokine gene IL-8 (CXCL-8) to distinguish temperature at 30mins but there are a small number of other NF-κB target genes such as NUAK2, NFKBIA, TNFAIP3 (A20) that could have been used instead (see [3]).

We use microarrays and RT–qPCR data to monitor the expression of these genes around the peak times of nuclear NF-κB at 0, 30, 130, 230 and 430 minutes (Fig 2a). We see that, if we know the expression levels of these genes at these times we can determine which of these multiple experiments was carried out (Fig 2b). Monitoring the gene expression at 30 minutes enables the identification of 5 distinct input signal combinations (unstimulated and the four combinations in the table) and also one additional one if observations at 130 minutes are included. This suggests that just from monitoring these three genes we obtain at least 2 bits of information.

thumbnail
Fig 2. Gene expression can identify different experimental conditions.

(a) (Left Column) The expression of the gene EGR1 in normal (37°C) and high (40°C) temperature and under continuous or pulsed TNFα treatment with pulses repeated every 100 minutes. The pulse length is 5 minutes except in the bottom row where the pulse length is indicated in the legend. (Middle Column) As Left Column except that the gene is COX-2 (PTGS2). (Right Column) As the first two rows of the Left Column except that the gene is IL-8 (CXCL-8). (b) A table showing which gene expression combinations identify which pairs of the input signal. The letters E, C and I indicate the genes EGR1, COX-2 and IL-8 respectively. The plus symbol (+) indicates high expression at 30 minutes and the minus symbol (-) indicates low expression at this time. Thus E+I indicates that EGR1 is highly expressed at 30 minutes and IL-8 is then at a low level, which implies that the system is pulsed with 5 minute pulses and the temperature is 37°C. The symbol +130 indicates high expression at 130 minutes.

https://doi.org/10.1371/journal.pcbi.1008076.g002

There are two substantial caveats to this observation. Firstly, a potential criticism is that the genes discussed might also be regulated by independent parallel pathways. However, we note that in Fig 2(b) we have restricted to very early observations at 30 minutes when the involvement of other pathways is unlikely. Secondly, we are using data from cell population assays such as microarrays rather than single cells where stochastic effects are important. However, the expression difference in the table are more than two logs so the overlap of the corresponding expression distribution should be small.

Despite these two caveats these results are highly suggestive and motivate a careful consideration of single cell multiplexing. Further information supporting these ideas is contained in Section 5 and Fig K in S1 Appendix.

Results

Decision-making and KL divergence

We now develop a mathematical approach that enables us to quantify multiplexing. We use this to show why current tightly coupled models of NF-κB cannot multiplex effectively and then explain how to modify these so that multiplexing is enabled.

Suppose we have s signals S1, …, Ss which in turn define the vector signal S = (S1, …, Ss). Consider a change in the signal from a base value S0 to S = S0 + δS where the change has size η = ∥δS∥. We ask whether, the response R has the capacity to distinguish which components of the signal have significantly changed i.e. to identify which components of δS are ≥ O(η). If it can we say the system can multiplex.

Mathematically the question of using the stochastic response R, which has probability distribution PS(R) = P(R|S), to distinguish input signals is related to hypothesis testing. If a change in input signal occurs (say from S0 to S = S0 + δS) and we wish to determine if the ith component was changed using only R, we need to be able to evaluate the hypothesis that R comes from PS rather than from a distribution of the form PS where S′ is any perturbation of S0 with the same ith component as S0. By the Neyman-Pearson lemma, the most powerful test of this hypothesis for a given false-positive error rate α is a test of the form λ(R) ≥ uα where is the log-likelihood ratio and the choice of α determines what threshold uα to use. The PS-mean of the log-likelihood ratio is by definition the Kullback–Leibler (KL) divergence, DKL(PSPS), of PS and PS distributions. The larger is the likelihood ratio, the more evidence we have in favour of signal S and against S′.

Multiplexing capacity

If DKL is too small then the most powerful test is expected to fail and hence other tests will not fair any better. Furthermore, as we wish to check whether the response R has the capacity to distinguish S from any signal S′ that has the i-th signal unchanged, we study how large is the where is the set of all such S′ signals. However, as S tends towards S0 thus decreasing the length l = l(S) = ∥SS0∥, this quantity decreases like l2, and therefore we scale it and define (1) The larger is, the easier it is to detect the change in the ith component.

To apply this so as to detect changes in any component of the signal S we consider (2) The larger this multiplexing capacity MX(S1, …Ss) is, the better the system at multiplexing the signals S1, …Ss (see also Sects. 2.1-2.4 in S1 Appendix).

Characterising multiplexing via the sensitivity matrix

While we cannot calculate this quantity in general, we can find an elegant solution in terms of the Fisher Information Matrix (FIM) when the changes in the signal are small, that is the third order terms and above are negligible. That is, we will calculate DKL(PSPS) up to terms that are O(max{∥SS′∥3, ∥SS03, ∥S′ − S03,}).

In our context, the FIM at S0 has entries where (S;R) = log P(R|S) is the log–likelihood function, ∂iℓ denotes the partial derivative with respect to the ith component Si and is the corresponding second derivative. These derivatives are evaluated at S0.

The FIM measures the sensitivity of PS(R) to a change δS in the signal S because, up to terms that are O(∥δS3) (see Sect. 2.1 in S1 Appendix),

Multiplexing sensitivity matrix.

One can associate to the FIM an s × s matrix that satisfies and certain optimality properties described in [17] (Sect. 2.5.2 in S1 Appendix). For j = 1, 2, …, s, we call the entries sij of the matrix s, the principal coefficients of sensitivity of the response R to the j-th signal Sj.

The multiplexing sensitivity matrix s describes the ability of the signalling system to multiplex at least locally in the following way. If sj, j = 1, …, s, are the columns of s, then

  1. denotes the linear subspace of spanned by the vectors , and
  2. n = n(i|i1, …, ik) denotes the component of si normal to the linear subspace i.e. si = u + n with u in and n orthogonal to . If i1, …, ik include all indices except j we use the notation n(i|ji) for n(i|i1, …, ik).

Firstly, up to third order terms (see Sect. 2.4 in S1 Appendix), and therefore the length of the normal component, n(i|ji), determines, at least locally, the capacity of the response R to distinguish the i-th from the rest of the considered signals. Secondly, there is an essentially unique reordering of the signal components as so that if vk = ∥n(ik|i1, …, ik−1)∥ then v1 ≥ ⋯ ≥ vs and the multiplexing capacities up to third order terms are (3) for all k = 1, …, s. All of these quantities can be rapidly calculated using the QR decompositions of submatrices of s made up from the relevant columns of s (see Sect. 2.4 in S1 Appendix).

This ordering of the set of signals provides a way to choose an optimal subset that can multiplex. That is, we can use the ordering i1, …, is and the associated multiplexing capacities , k = 1, …, s, to identify the subset of signals with the largest number of elements, k, that has multiplexing capacity , for m an appropriate threshold (e.g. the minimum DKL level for the change to be detectable in a given system of interest, see also Sect. 2.4.1 in S1 Appendix).

Multiplexing capacities of a model

In regulatory and signalling systems, the values of two parameters, say and , may differ by an order of magnitude or more. Therefore, when discussing sensitivities it is usually not appropriate to consider the absolute changes in the parameters , but instead to consider the relative changes. A good way to do this is to introduce new parameters because absolute changes in θj correspond to relative changes in . Then, for small changes to the parameters, and so the changes δθj are scaled and non-dimensional. When discussing the multiplexing capacities of a model below we will always use these scaled parameters θj and henceforth when we refer to parameters θi we mean these scaled parameters and we drop the tilde.

As mentioned above the typical situation is where the signals Si change the parameters θj, j = 1, …, s, of the model so that θ = θ(S). In this case we show in Sect. 2.5.3 in S1 Appendix how one can relate the multiplexing capacity of the signals to the multiplexing capacity of the parameters. For the latter we regard the parameters as signals and calculate their multiplexing capacities where the parameters θi have been reordered so that for all m < s, . We call these the multiplexing capacities of the model.

Knowing the multiplexing capacities of a model is important because, using equation (15) of Sect. 2.5.3 in S1 Appendix, we can tell how well the system can multiplex any signals changing these parameters. In particular,

  1. if the MXm very rapidly decrease with m then the system is not able to multiplex through these parameters; and
  2. if is large and there are signals S1, …, Sr which change the parameters in that , then these signals will have good multiplexing properties provided the matrix d θ(S) (which is the derivative of θ with respect to S evaluated at the base value of S) is well conditioned (e.g. if detd θ(S) is not too small or big).

Tightly coupled models of NF-κB cannot multiplex effectively

One might expect that a dynamical system with many parameters, such as NF-κB (see Fig 3(a)), would have the flexibility to multiplex effectively. However, it has been observed that for a large class of deterministic models of regulatory and signalling systems of the sort that we are considering, the deterministic analogue of the FIM for the model parameters has rapidly decreasing eigenvalues [8, 3135]. A similar result was shown for stochastic models of the circadian clock in [17, 36]. This implies that the effects of changing different parameters are highly correlated making it hard to recognise which parameter was changed.

thumbnail
Fig 3. Diagrams of the base NF-κB model in [48] and its modification mNF-κB model.

(a) The main reactions following a TNFα stimulus according to the model in [48] (base model); (b) The main reactions of the mNF-κB model following a TNFα stimulus, and when the modification signal S2 is constantly transmitted. For the m2NF-κB model, the S2 signal controls NF-κB modification jointly (AND logic) with the TNFα signal through the active IKK molecules (dashed line).

https://doi.org/10.1371/journal.pcbi.1008076.g003

We see in Fig 4(a) that such a rapid decline in the eigenvalues of the FIM is the case for the base model considered here. They decay with an exponential rate and the second singular value is already less than 1% of the first one.

thumbnail
Fig 4. Comparisons of the multiplexing capacities and sensitivities of the base NF-κB and mNF-κB models.

(a) The singular values σi of the FIM for the base model and the much larger singular values for the mNF-κB model. (b) The multiplexing capacities of the base and mNF-κB models. The parameters with the largest multiplexing capacities correspond to the scaled version of the parameters TNFKB (total amount of NF-κB molecules) and TNFα dose for the base model and TNFKB, pd1 (reverse modification rate of NF-κB) and S2 (signal; treated as parameter here) for the mNF-κB model. (c) The principal sensitivity coefficients of the mNF-κB model. Larger values indicate higher sensitivity of the mNF-κB model to changes in the value of the corresponding parameter. (d,e) Realisations (n = 1000) of the pcLNA distributions of the base and mNF-κB model respectively at three times chosen to correspond to the first three peaks of the deterministic limit (Ω → ∞) of the model. In each case the 6 clusters correspond to different scaled parameter values which are either the scaled base value θ0 or the scaled parameter vector θ0 + δVj, j = 1, …, 5 (δ = 0.1) where V1, …, V5 are the eigendirections of the FIM corresponding to the 5 largest singular values of each model. Thus each cluster of points corresponds to a different principal component. Notice how much better the clusters are separated in the mNF-κB model compared to the base NF-κB model. For the mNF-κB model the base and first three principal component perturbations (black, red, green and blue) are effectively completely separated.

https://doi.org/10.1371/journal.pcbi.1008076.g004

But and therefore the eigenvalues of the FIM are the squares of the singular values σi of s (Sect. 2.3 in S1 Appendix). Consequently the singular values σi of s rapidly decrease and, since v1vkσ1σk for all ks with equality when k = s (Theorem 3.3.2 of [37], see Sect. 1 in S1 Appendix) the same is true for the multiplexing capacities , j = 1, …, k. Using Eq (3) we see that the number of signals that can multiplex well must be very small. Fig 4(b) shows how fast the multiplexing capacities decrease for our base model and identifies the parameters through which signals can multiplex more effectively i.e. these are the parameters that the signals should move if the signals are to be effective.

Increased multiplexing via additional regulated NF-κB states

Regulation of the NF-κB pathway is enabled by multiple post-translational modifications that control the activity of the core components of NF-κB signaling. In particular, the RelA NF-κB subunit undergoes reversible modifications such as phosphorylation, ubiquitination, and acetylation that can affect its transcriptional functions [3844]. Indeed many modification sites in RelA have been identified as having either an enhancing, inhibitory, or modulatory effect on NF-κB transcriptional activity in a gene-specific manner [3841, 43]. A further potentially regulated step that could differentially control individual gene expression is the hetero- and homo-dimerisation of the NF-κB Rel proteins [45, 46]. Therefore, in considering the nature of biological mechanisms that could underlie multiplexing of information by the NF-κB system, it is natural to consider modifications that create additional regulated NF-κB states that can affect the transcription of NF-κB target genes.

We consider one of the simplest modifications of the base model that can enable more effective multiplexing. In this modified model, which we call mNF-κB, the cytoplasmic NF-κB is reversibly modified by an input signal S2 that is independent to the TNFα signal. For example, S2 might be an environmental signal such as temperature or pressure that strongly affects the activity of a kinase or other molecular processes. Such temperature effects have been studied extensively for the circadian clock in the context of temperature compensation [47] and more recently for the NF-κB system [3]. Moreover, in Fig 2 we see substantial temperature effects on the genes considered there.

The modified form of NF-κB, mNF-κB, competes with the unmodified form for binding of IκBα but otherwise is subject to the same reactions (see Fig 3(b) and Sect. 3.3 in S1 Appendix). Importantly, mNF-κB can activate, inhibit, or modulate the transcription of target genes and their differential expression can potentially reveal the levels of the S2 signal. The mathematical analysis of the stochastic version of the mNF-κB model confirms this in the following ways.

Firstly, the singular values of the FIM are overall increased, and, importantly, there are now two large singular values rather than one (Fig 4(a)). Secondly, the multiplexing capacities are increased substantially in the mNF-κB model with three of them significantly above the second for the NF-κB model. As explained above in the section “Multiplexing capacities of a model” and Sect. 2.5.1 in S1 Appendix this means that this model supports multiplexing of three signals through these three parameters. That the multiplexing capacities are large for the parameters related to the modification confirms that the extra sensitivity arises from the addition of this modification (see Fig 4(d)).

Note that the results presented in Fig 4 are derived for the probability distributions of stochastic trajectories of the system observed at 9 timepoints (see Sect. 4.9 in S1 Appendix). If instead only two time-points are considered, that is 10 mins before and at the expected time of the first peak of nuclear NF-κB concentration, the base NF-κB model is not largely affected, but the mNF-κB presents a clearly less prominent increase of the singular values (see Fig H in S1 Appendix). This suggests that while the dynamical behaviour of the system does not in itself enable higher multiplexing capacity, it can greatly enhance multiplexing in a system that has the ability to multiplex.

The greater sensitivity of the mNF-κB model compared to the base model is also reflected in Fig 4(c). We see that the nuclear concentrations of NF-κB are much more affected in the mNF-κB model by changes in the signal. This clearly provides much greater ability in modulating gene expression according to different signals (see next section).

In this example we used one of the most simple and generic modifications where an external environmental signal such as temperature affects an internal parameter. In Sect. 3.5 in S1 Appendix we also consider another modification where the variation in the internal parameter is caused by a noisy pathway. Clearly, in this case the multiplexing capacity depends on the mechanism of this modification pathway and its information carrying effectiveness. Despite the increased noise levels, the multiplexing capacity of this alternative model can also be significantly larger than the base model.

Reproducing multiplexing in the EGR1-COX-2 example

To further illustrate multiplexing, we now consider how to modify the signalling system so as to be able to reproduce the multiplexing behaviour seen in the EGR1-COX-2 gene expression data in Fig 2. The same principle can be extended to reproduce the expression of the 8 genes presented in Fig K in S1 Appendix, but this is beyond our scope, and more data will be necessary for validation. Furthermore, we are not claiming that this is the true underlying biological mechanism but are using this example to illustrate how the NF-κB signalling system can multiplex different signals through gene regulation. This is clearly not possible under the structural constraints of the base model because: (a) the base NF-κB model reacts to pulses of TNFα stimuli by nearly identical (forced) oscillations and therefore it cannot explain the difference between the early and late expression of EGR1 and COX-2, and (b) the differences in the base model between the response to short and long pulse are extremely small and can hardly explain the differences in EGR1 early response between the different pulse lengths.

The system is modified to include a reversible modification of NF-κB molecules in the cytoplasm. The NF-κB modification is jointly promoted by the TNFα stimulus through the IKK module and the independent signal S2 (see Fig 3(b)). Pulses of TNFα cause bursts of NF-κB nuclear translocations, but also higher levels of the modified NF-κB. The reverse modification is independent of S2 and TNFα. Apart from TNFα promoting the NF-κB modification, this model which we call m2NF-κB is the same as the mNF-κB model (see Fig 3(b) and Sect. 3 in S1 Appendix).

The m2NF-κB model postulates that NF-κB activates the transcription of EGR1, which is inhibited by the mNF-κB, while the reverse regulation is imposed on COX-2. Using our approach to stochastic simulation outlined next, we can calculate the confidence limits for COX-2 and EGR1 under the various pulsing protocols (see Fig 5) using n = 1000 trajectories simulated as described in the next section (see also Sect. 4.6 in S1 Appendix). Fig 5(a) provides the mean time-trajectories (and 10 samples) at the same times observed using microarray and qPCR in Fig 2. The introduction of the additional regulatory states of NF-κB allows us to reproduce the experimentally observed profile.

thumbnail
Fig 5. Stochastic simulation of the mNF-κB model II describing the regulation of EGR1 and COX-2 genes for the same types of TNFα stimulation as in Fig 2.

(a) A sample of n = 10 realisations (dots) and the sample mean (straight line) of simulated trajectories under the different TNFα stimuli (see legend) at the same times (t = 0, 30, 130, 230, 430min) as observations in Fig 2; (b) 95% confidence envelopes of the copy number of (unbound) nuclear NF-κB and mNF-κB molecules, and EGR1, COX-2 mRNA copies under the different TNFα stimuli derived using stochastic simulations (n = 1000) under the different TNFα stimuli (see legend). The base model cannot reproduce the observed sensitivity to the different pulse lengths.

https://doi.org/10.1371/journal.pcbi.1008076.g005

Stochastic dynamics of NF-κB

The base model used in our analysis is a stochastic reaction network that describes the oscillatory response of the NF-κB system under stimulation by TNFα. It is a slight modification of the system model in [48]. In our version of the model, after adjustments to the rate equations, concentrations are all expressed in terms of the same volume Ω, taken to be Avogadro’s number in the appropriate molar units multiplied by the volume of the cell in appropriate units so that Ω has units L/nM (Sect. 3.3 in S1 Appendix). The original model is written in terms of nuclear and cytoplasmic concentrations. Clearly, it is straightforward to convert between the two models (see Sect. 3.3 in S1 Appendix).

We use the pcLNA stochastic version of this model [17] that allows us to derive analytical expressions for the FIM and system sensitivity matrix s and to rapidly simulate the system with high accuracy (see Fig 6 and Sect. 4.2 in S1 Appendix). The stochastic model considered here converges to the published deterministic model of [48] as Ω → ∞. We believe that the ability of our method to calculate important information-theoretic quantities such as these for a large fully stochastic model is a significant new development in itself.

thumbnail
Fig 6. The pcLNA stochastic model and the channel capacity of the NF-κB model.

(a) The pcLNA model uses the stability of the probability distributions of stochastic oscillatory systems on the transversal sections, Sx, of a given phase, x, of the system’s deterministic solution. (b) The pcLNA probability distributions on those transversal sections match very well the empirical distributions derived by SSA. Here the comparison is done using the Kolmogorov-Smirnov (KS) test at the first 4 peaks of NF-κB model. The corresponding histograms for two of the largest observed KS values are also displayed to illustrate the nearly perfect match of the two distributions even in the case of the largest KS distances recorded here. (c) The pcLNA simulations also match very well the SSA simulations of the NF-κB model which is much slower (see CPU (average) time for a single simulation). (d) Estimation of the channel capacity using the pcLNA simulation algorithm with added noise on the Total number of NF-κB molecules.

https://doi.org/10.1371/journal.pcbi.1008076.g006

The mNF-κB model that includes NF-κB modification is also simulated and analysed using pcLNA (see Sect. 3.4 in S1 Appendix). For the simulation of downstream genes that are regulated by NF-κB (see next section) we use the Stochastic Simulation Algorithm (SSA) [49]. This is because the relevant distribution for the gene expression is far from being Gaussian and therefore it is not appropriate to apply the pcLNA directly to this subsystem. Since this part of the system involves relatively few molecules the combined system can be simulated rapidly. The SSA is also used for comparisons to pcLNA in Fig 6. In Fig 6, we show that pcLNA accurately approximates the SSA simulations of the model in [48], which was calibrated to experimental data.

Capacity of scalar channels

The results above raise the question of whether our models are compatible with the channel capacity seen in previous publications. We can use the base model to compare its behaviour with that discussed in [19, 21]. In these papers there was no attempt to control cell size or consideration of the total amount of NF-κB (see also [23] which discusses such issues). We therefore allow these quantities to vary with the variation being drawn from a log-Normal distribution as described in Sect. 4.4 in S1 Appendix.

We study the case where S is the level of the continuous TNFα stimulation (the parameter dose) and the response R is the level of nuclear NF-κB at q different phases including its first peaks and troughs. Fig 6(d)(i) shows the estimated capacity as a function of q. We also estimate the channel capacity for response R the nuclear concentration at t = 30min after initiating continuous TNFα stimulation (Fig 6(d)(ii)).

The model reproduces the rather limited channel capacity seen in [19, 21] with estimated carrying capacities in the region of one bit. The exact value is not important because this is subject to our estimates of Ω, the total concentration of NF-κB molecules and other parameters derived in [48]. A similar result can be obtained using the model in [30] (see Sect. 4.5 in S1 Appendix). It is worth emphasizing that the limited channel capacity is observed for a scalar signal and scalar response.

Discussion

Cells present a very different context from that of traditional communications channels. The genetic and epigenetic information contained in the genome is translated by molecular interactions into dynamical processes. Described by dynamical interaction networks, these stochastic dynamical processes effectively move information from one system to another by regulating the probability distributions of their component molecules. Therefore, it is unclear whether the classical tools are always the most appropriate and it is likely that a much more extensive information toolbox is needed. New ideas about stochasticity and information are needed to understand how cells respond to dynamic environments so as to ensure appropriate cellular responses with high probability when they are using biochemistry that itself is very noisy.

Using such information theoretic tools we suggest a new insight into the way in which signalling systems transmit information. We mentioned above that recent research [19, 21] has shown that the channel TNFα level → nuclear NFKB abundance has a relatively low channel capacity. This raises the question of how our results fit with this. To some extent this is answered by the results in the section entitled “Capacity of scalar channels” where we show that our systems are tuned so as to reproduce this. Clearly, if we ignore noise and use a deterministic system we can make any such channel have as large a capacity as we want so it is important in our work to use reasonable levels of stochasticity.

We suggest that there is coherent picture emerging here where although signalling systems may be rather limited in the way that they transmit any scalar signal (as above for TNFα level), they are well designed to transmit multi-dimensional signals. There are two main reasons why when considering information flows in signalling one wishes to consider gene responses that are multidimensional. The first is that transmitting a signal via multiple receivers enables one to reduce the effects of noise. The second which is of central concern here is that it enables complex non-binary decisions. However, to make use of all these dimensions it is necessary that the input signal S has multiple dimensions because otherwise, if R is d-dimensional, the mean of P(R|S) is constrained to a 1d curve in d-dimensional space. This would mean that to obtain multiplexing or higher channel capacity one would have to use changes in the variance of P(R|S) with S to detect changes which seems very unlikely to be effective. Indeed, to use all d dimensions one needs dim Sd.

We envisage that in this multidimensional situation it may well be the case that the scalar channels SiR each have very low capacity as is the case in [19, 21] but that the full system SR is able to multiplex so as to enable complex decisions and has a significantly higher capacity. Thus, by using multi-dimensionality the system can use multiple low-capacity components to produce a high capacity system.

A related issue concerns the role of dynamics in information transfer including the suggestion that dynamic systems such as oscillating ones can transmit greater amounts of information compared to static/equilibrium systems [48, 5055]. Our examples, also suggest why an oscillating system can use multiplexing to transmit more information than equilibrium systems. In these we see that signals that affect protein modification states or other aspects such as dimerization or binding partners can be good for multiplexing. In an equilibrium system the probability distribution describing how these states are distributed will be stationary in time. On the other hand in an oscillatory system these states can have a non-trivial temporal structure (e.g. oscillating) as catalysts of modifications can be activated and deactivated by interaction with the oscillations. This suggests a clear advantage for oscillating systems for information transfer.

Methods

Cell culture and reagents

Experiments were performed on human neuroblastoma SK-N-AS cells cultured in Modified Eagles Medium supplemented with 1% non-essential amino acids (Sigma-Aldrich) and 10% foetal bovine serum (Gibco). Cells were maintained in a humidified 5% CO2 incubator at 37°C.

Microarray and RT-qPCR experiments

SK-N-AS cells were plated at a density of 500,000 cells per dish. Time course experiments were carried out 24 hours later. Cells were transferred to the three temperatures, 34°C, 37°C and 40°C for 1h and then (at time 0) stimulated with 10ng/ml of TNFα. TNFα was added continuously, for a single pulse of varying durations (2.5, 5, 10, 20, 40 minutes) or as repeated 5 minute pulses of 60/100/200 minutes intervals. Diclofenac (300, 500 μg/ml; Sigma Aldrich) was added 1.5h prior to TNFα. Measurements were taken at 0, 15, 30, 130, 230 and 430 minutes. Cells were lysed, total RNA was extracted and RT-qPCR and microarray experiments were performed as described previously [3]. Primer sequences used were COX-2 left—gcaataacgtgaagggctgt, right—cgggaagaacttgcattgat, EGR1 left—ttcccttcctcagctgtcac, right—tgtcctgggagaaaaggttg. Data on the effect of temperature on TNFα-induced gene expression generated previously were used in this study [3].

Supporting information

S1 Appendix. Supporting results, details of mathematical analysis, and description of the computational algorithms and the models used in the main manuscript.

https://doi.org/10.1371/journal.pcbi.1008076.s001

(PDF)

References

  1. 1. Hao S, Baltimore D. The stability of mRNA influences the temporal order of the induction of genes encoding inflammatory molecules. Nat Immunol. 2009; 10(3): 281–8. pmid:19198593
  2. 2. Zhang Q, Lenardo MJ, Baltimore D. 30 Years of NF-κB: A Blossoming of Relevance to Human Pathobiology. Cell. 2017; 168(1-2) 37–57. pmid:28086098
  3. 3. Harper CV, Woodcock DJ, Lam C, Garcia-Albornoz M, Adamson A., Ashall L, et al. Temperature regulates NF-κB dynamics and function through timing of A20 transcription. PNAS. 2018; 115(22) E5243–E5249. pmid:29760065
  4. 4. Cover TM, Thomas JA. Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing). New York, NY, USA: Wiley-Interscience; 2006.
  5. 5. Cacuci DG. Sensitivity and Uncertainty Analysis: Theory. I. Boca Raton, Florida: Chapman & Hall; 2003.
  6. 6. Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S. Global Sensitivity Analysis: The Primer. 2008; John Wiley & Sons.
  7. 7. Saltelli A, Tarantola S, Campolongo F, Ratto M. Sensitivity analysis in practice: a guide to assessing scientific models (Vol. 1). New York: Wiley; 2004.
  8. 8. 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), 1871–8. pmid:17922568
  9. 9. Hwang JT, Dougherty EP, Rabitz S, Rabitz H. The Green’s function method of sensitivity analysis in chemical kinetics. The Journal of Chemical Physics. 1978; 69(11) 5180–5191.
  10. 10. Ljung L. Perspectives on system identification. Annual Reviews in Control 34.1. 2010: 1–12.
  11. 11. Plyasunov S, Arkin AP, Efficient stochastic sensitivity analysis of discrete event systems. J. Comput. Phys. 2007; 221, 724–738.
  12. 12. Wolf ES, Anderson DF, A finite difference method for estimating second order parameter sensitivities of discrete stochastic chemical reaction networks. J. Chem. Phys. 2012; 137, 224112.
  13. 13. Gupta A, Rathinam M, Khammash M, Estimation of parameter sensitivities for stochastic reaction networks using tau-leap simulations, SIAM J. Numer. Anal. 2014; 56, 1134–1167.
  14. 14. Komorowski M, Costa MJ, Rand DA, Stumpf MPH, Sensitivity, robustness, and identifiability in stochastic chemical kinetics models, PNAS. 2011; 108(21) 8645–8650. pmid:21551095
  15. 15. Villaverde AF, Banga JR, Dynamical compensation and structural identifiability of biological models: Analysis, implications, and reconciliation. PLoS Computational Biology. 2017; 13(11): e1005878. pmid:29186132
  16. 16. Sontag ED. Dynamic compensation, parameter identifiability, and equivariances. PLoS Computational Biology. 2017; 13(4): e1005447. pmid:28384175
  17. 17. Minas G, Rand DA. Long-time analytic approximation of large stochastic oscillators: Simulation, analysis and inference PLoS Comput Biol. 2017; 13(7):e1005676. pmid:28742083
  18. 18. Bowsher CG, Voliotis M, Swain PS. The Fidelity of Dynamic Signaling by Noisy Biomolecular Networks. PLoS Computational Biology. 2013; 9(3):e1002965. pmid:23555208
  19. 19. Cheong R, Rhee A, Wang C, Nemenman I, Levchenko A. Information Transduction Capacity of Noisy Biochemical Signaling Networks. Science. 2011; 334(6054):354–8. pmid:21921160
  20. 20. Perkins TJ, Swain PS. Strategies for cellular decision-making. Mol Syst Biol. 2009; 5:326. pmid:19920811
  21. 21. Selimkhanov J, Taylor B, Yao J, Pilko A, Albeck J, Hoffmann A, Tsimring L, Wollman R. Systems biology. Accurate information transmission through dynamic biochemical signaling networks. Science. 2014; 346(6215):1370–3. pmid:25504722
  22. 22. Voliotis M, Perrett RM, McWilliams C, McArdle CA, Bowsher CG Information transfer by leaky, heterogeneous, protein kinase signaling systems. PNAS. 2014; 111(3):E326–33. pmid:24395805
  23. 23. Zhang Q, Gupta S, Schipper DL, Kowalczyk GJ, Mancini AE, Faeder JR, Lee REC. NF-κB Dynamics Discriminate between TNF Doses in Single Cells. Cell Systems. 2017; 5(6), 638–645.e5. pmid:29128333
  24. 24. Kellogg RA, Tay S. Noise facilitates transcriptional control under dynamic inputs. Cell. 2015; 160(3), 381–392. pmid:25635454
  25. 25. Kellogg RA, Tian C, Etzrodt M,Tay S, Ryan Kellogg AA, Tian C, Tay S. Cellular Decision Making by Non-Integrative Processing of TLR Inputs. Cell Reports. 2017; 19(1), 125–135. pmid:28380352
  26. 26. Lane K, Van Valen D, DeFelice MM, Macklin DN, Kudo T, Jaimovich A, Covert MW. Measuring Signaling and RNA-Seq in the Same Cell Links Gene Expression to Dynamic Patterns of NF-κB Activation. Cell Systems. 2017; 4(4), 458–469.e5. pmid:28396000
  27. 27. Lee REC, Qasaimeh MA, Xia X, Juncker D, Gaudet S. NF-κB signalling and cell fate decisions in response to a short pulse of tumour necrosis factor. Scientific Reports. 2016; 6, 39519.
  28. 28. Tudelska K, Markiewicz J, Kochaczyk M, Czerkies M, Prus W, Korwek Z, Lipniacki T. Information processing in the NF-κB pathway. Scientific Reports. 2017; 7(1), 15926. pmid:29162874
  29. 29. Suderman R, Bachman JA, Smith A, Sorger PK, Deeds EJ. Fundamental trade-offs between information flow in single cells and cellular populations. PNAS. 2017; 114(22), 5755–5760. pmid:28500273
  30. 30. Tay S, Hughey J, Lee T, Lipniacki T, Quake S, Covert M. Single-cell NF-B dynamics reveal digital activation and analogue information processing. Nature. 2010; 7303(466):267–271.
  31. 31. Rand DA, Shulgin BV, Salazar D, Millar AJ. Design principles underlying circadian clocks. J R Soc Interface. 2004; 1(1) 119–30. pmid:16849158
  32. 32. Rand DA, Shulgin BV, Salazar D, Millar AJ. Uncovering the design principles of circadian clocks: Mathematical analysis of flexibility and evolutionary goals. J. Theor. Biol. 2006; 238(3):616–35. pmid:16111710
  33. 33. Rand DA. Mapping the global sensitivity of cellular network dynamics: Sensitivity heat maps and a global summation law. J. R. Soc. Interface. 2008; 6(5) Suppl 1:S59–69.
  34. 34. Waterfall JJ, Casey FP, Gutenkunst RN, Brown KS, Myers CR, Brouwer PW, Elser V, Sethna JP. Sloppy-Model Universality Class and the Vandermonde Matrix Phys. Rev. Lett. 2006; 97(15):150601. pmid:17155311
  35. 35. Transtrum MK, Machta BB, Sethna JP. Why are Nonlinear Fits to Data so Challenging? Phys. Rev. Lett. 2010; 104(6):060201. pmid:20366807
  36. 36. Minas G, Rand DA. Parameter sensitivity analysis for biochemical reaction networks. Mathematical Biosciences and Engineering. 2019; 16(5): 3965–3987. pmid:31499645
  37. 37. Horn RA, Johnson CR Topics in Matrix Analysis. Cambridge: Cambridge University Press; 1991.
  38. 38. Perkins ND. Post-translational modifications regulating the activity and function of the nuclear factor kappa B pathway. Oncogene. 2006; 25(51):6717–6730. pmid:17072324
  39. 39. Christian F, Smith E, Carmody R. The Regulation of NF-κB Subunits by Phosphorylation. Cells. 2015; 5(1):12, mar 2016.
  40. 40. Collins PE, Mitxitorena I, Carmody RJ. The Ubiquitination of NF-κB Subunits in the Control of Transcription. Cells. 2016; 5(2):23.
  41. 41. Bo Huang, Xiao-Dong Yang, Acacia Lamb, and Lin-Feng Chen. Posttranslational modifications of NF-kappaB: another layer of regulation for NF-kappaB signaling pathway. Cellular signalling. 2019; 22(9):1282–90.
  42. 42. Lu X, Yarbrough WG. Negative regulation of RelA phosphorylation: Emerging players and their roles in cancer. Cytokine & Growth Factor Reviews. 2015; 26(1):7–13.
  43. 43. Wang B, Wei H, Prabhu L, Zhao W, Martin M, Hartley AV, Tao L. Role of Novel Serine 316 Phosphorylation of the p65 Subunit of NF-κB in Differential Gene Regulation. The Journal of biological chemistry. 2015; 290(33):20336–47. pmid:26082493
  44. 44. Zhang L, Shao L, Creighton CJ, Zhang Y, Xin L, Ittmann M, Wang J. Function of phosphorylation of NF-kB p65 ser536 in prostate cancer oncogenesis. Oncotarget. 2015; 6(8):6281–94. pmid:25749044
  45. 45. Tsui R, Kearns JD, Lynch C, Vu D, Ngo KA, Basak S, Ghosh G, Hoffmann A. IκBβ enhances the generation of the low-affinity NFκB/RelA homodimer. Nat Commun. 2015; 6: 7068. pmid:25946967
  46. 46. Bosisio D, Marazzi I, Agresti A, Shimizu N, Bianchi ME, Natoli G. A hyper-dynamic equilibrium between promoter-bound and nucleoplasmic dimers controls NF-κB-dependent gene activity. The EMBO Journal. 2006; 25(4):798–810. pmid:16467852
  47. 47. Hastings JW. & Sweeney BM On the mechanism of temperature independence in a biological clock. PNAS. 1957; 43:804–811. pmid:16590089
  48. 48. Ashall L, Horton CA, Nelson DE, Paszek P, Harper CV, Sillitoe K, et al. Pulsatile stimulation determines timing and specificity of NF-κB-dependent, transcription. Science 2009; 324(5924):242–6. pmid:19359585
  49. 49. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 1977; 81(25):2340–61.
  50. 50. Cai L, Dalal CK, Elowitz MB. Frequency-modulated nuclear localization bursts coordinate gene regulation. Nature. 2008; 455(7212): 485–90. pmid:18818649
  51. 51. Lahav G, Rosenfeld N, Sigal A, Geva-Zatorsky N, Levine AJ, Elowitz MB, Alon U. Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat Genet. 2004; 36(2):147–50. pmid:14730303
  52. 52. Nelson DE, Ihekwaba AE, Elliott M, Johnson JR, Gibney CA, Foreman BE, et al. Oscillations in NF-κB signaling control the dynamics of gene expression Science. 2004; 306(5696):704–8. pmid:15499023
  53. 53. Shankaran H, Ippolito DL, Chrisler WB, Resat H, Bollinger N, Opresko LK, Wiley HS Rapid and sustained nuclear-cytoplasmic ERK oscillations induced by epidermal growth factor. Mol. Syst. Biol. 2009; 5:332. pmid:19953086
  54. 54. Yoshiura S, Ohtsuka T, Takenaka Y, Nagahara H, Yoshikawa K, Kageyama R. Ultradian oscillations of Stat., Smad., and Hes1 expression in response to serum. PNAS. 2007; 104(27):11292–7. pmid:17592117
  55. 55. Purvis JE, Lahav G. Encoding and decoding cellular information through signaling dynamics. Cell. 2013; 152(5):945–956. pmid:23452846