1. Computational and Systems Biology
  2. Genomics and Evolutionary Biology
Download icon

NF-κB oscillations translate into functionally related patterns of gene expression

  1. Samuel Zambrano
  2. Ilario De Toma
  3. Arianna Piffer
  4. Marco E Bianchi Is a corresponding author
  5. Alessandra Agresti Is a corresponding author
  1. San Raffaele Scientific Institute, Italy
  2. San Raffaele University, Italy
Research Article
Cited
15
Views
2,377
Comments
0
Cite as: eLife 2016;5:e09100 doi: 10.7554/eLife.09100

Abstract

Several transcription factors (TFs) oscillate, periodically relocating between the cytoplasm and the nucleus. NF-κB, which plays key roles in inflammation and cancer, displays oscillations whose biological advantage remains unclear. Recent work indicated that NF-κB displays sustained oscillations that can be entrained, that is, reach a persistent synchronized state through small periodic perturbations. We show here that for our GFP-p65 knock-in cells NF-κB behaves as a damped oscillator able to synchronize to a variety of periodic external perturbations with no memory. We imposed synchronous dynamics to prove that transcription of NF-κB-controlled genes also oscillates, but mature transcript levels follow three distinct patterns. Two sets of transcripts accumulate fast or slowly, respectively. Another set, comprising chemokine and chemokine receptor mRNAs, oscillates and resets at each new stimulus, with no memory of the past. We propose that TF oscillatory dynamics is a means of segmenting time to provide renewing opportunity windows for decision.

https://doi.org/10.7554/eLife.09100.001

eLife digest

The process of producing useful biological molecules from genes – known as gene expression – is not always simple. Many genes are part of complex circuits, some of which show regular patterns of activity in response to an environmental cue. For example, the expression of some genes is tied to the 24-hour daily cycle of light and dark.

Transcription factors are proteins that control gene activation and expression, and some transcription factors periodically move in and out of the cell’s nucleus – the compartment of an animal cell that houses the vast majority of the genetic material. This behavior is known as oscillation. A transcription factor called NF-κB oscillates, changing between an inactive form outside of the nucleus and an active form inside. NF-κB plays important roles in inflammation and cancer, and is activated by cues from outside the cell. Some of the genes that the active form of NF-κB activates then produce molecules that inactivate NF-κB, thus helping to establish the oscillations.

The benefits of the oscillations are not clear. However, recent studies suggest that environmental cues can cause small perturbations that gradually adjust the rate at which the oscillations occur, and in doing so, synchronize the oscillations amongst neighboring cells.

By using embryo cells from genetically engineered mice, Zambrano et al. investigated how NF-κB oscillations get synchronized. The experiments showed that the activity of the NF-κB protein and the expression of the genes it controls synchronize across neighboring cells whenever the external environmental perturbations come in pulses. However, once the pulsed cues stop, this synchronization is quickly lost. In essence, the cells reset after each environmental cue with no memory of previous episodes of NF-κB activity.

Further work revealed that the expression of the genes controlled by NF-κB also cycles and resets with each new environmental cue. However, the products of these genes accumulate in three different ways. Some accumulate quickly; some accumulate at a slow and steady pace; and some oscillate in amount, and this amount resets once the environmental cue has stopped. Each of these classes of gene products can be related to specific cell behaviors that activate sequentially on well-defined time schedules.

Overall, Zambrano et al. suggest that the ability of NF-κB to reset its activity with each new environmental cue gives cells the opportunity to pause and adjust course.

Zambrano et al. now plan to explore what happens to NF-κB synchronization in different cell types exposed to a collection of inflammatory stimuli. Along the same line, it will be worth exploring NF-κB behavior in cancer cells, where NF-κB activity is often out of control and drives unrestrained cell proliferation. These studies would contribute to a deeper understanding of cancer biology and to the identification of new treatments for the disease.

https://doi.org/10.7554/eLife.09100.002

Introduction

Genetic circuits are instrumental for cells to provide adequate transcriptional responses to different external and internal stimuli, but we are still far from a complete understanding of how they work. The growing availability of single-cell measurements is bringing unprecedented insights into the dynamics of transcriptional responses: although it was traditionally thought that gene circuits should provide a temporally stable response to constant external stimuli, there is increasing evidence that the dynamics of gene circuits is more elaborate, and is often characterized by pulses of activity (Levine et al., 2013). A special case is oscillatory dynamics, in which the activity of some element in the genetic circuit varies periodically, and so does the output of the circuit. Paradigmatic oscillatory dynamics are associated with circadian clocks, which are present in a wide variety of organisms that synchronize their activities to the periodic variation in external daylight (Bell-Pedersen et al., 2005). However, oscillatory dynamics with periods far shorter than 24 hr have been reported in genetic circuits of bacteria (Suel et al., 2007), yeast (Cai et al., 2008; Hao and O'Shea, 2012) and mammalian cells (Geva-Zatorsky et al., 2006; Hoffmann et al., 2002; Larson et al., 2013; Nelson et al., 2004; Shankaran et al., 2009).

The dynamics of transcription factors in gene circuits can play a key role in information transmission through biochemical networks (Selimkhanov et al., 2014) by selectively modulating the expression of different genes (Purvis and Lahav, 2013). A paradigmatic example is how different p53 dynamics can selectively activate transcriptional programs that commit the cell to different cell fate decisions (Purvis et al., 2012).

A property of free oscillating systems (or free oscillators, as analogy to simple mechanical oscillators such as the pendulum) is entrainment, by which an oscillator can gradually modulate its phase and frequency thanks to a small perturbation (forcing) that is itself periodic and oscillating with a period resonant with the intrinsic period of the oscillator (Pikovsky et al., 2003). A common forcing can also lead to the synchronous dynamics of multiple oscillators and thus to collective dynamical states. Entrainment was successfully reproduced in synthetic genetic oscillators, in which single bacterial cells (each expressing a fluorescent reporter protein in response to a synthetic genetic circuit) were entrained collectively to an oscillating provision of arabinose (Mondragon-Palomino et al., 2011). The forcing to one oscillator can also be provided by other oscillators, and this coupling can lead to the emergence of different collective dynamical states characterized by different synchronous dynamics (Pikovsky et al., 2003). Inter-cellular coupling has indeed been exploited to genetically engineer synchronous quorum sensing genetic oscillators in bacteria (Danino et al., 2010). On the other hand, intra-cellular coupling leads to locked oscillatory states for different cell oscillators, as recently shown for the circadian rhythm and the cell cycle (Bieler et al., 2014).

In mammalian cells, NF-κB is a typical transcription factor that displays intrinsic oscillatory behavior. NF-κB plays key roles in inflammation, immune responses, development and cancer (Chaturvedi et al., 2011; Ghosh and Hayden, 2008; Karin, 2006; Ledoux and Perkins, 2014; Naugler and Karin, 2008). Strictly speaking, NF-κB is a family of dimers encoded by 5 different genes, but in what follows we will refer to p65/RELA independently of the dimer it forms. In resting cells, p65 exists mostly within a cytoplasmic complex bound to the IκB inhibitors (Hoffmann et al., 2002). Inflammatory signals like tumor necrosis factor alpha (TNF-α) or lipopolysaccharide (LPS) induce phosphorylation of IκB proteins by IKKs–upstream kinases in the signaling pathway–, ubiquitination and degradation of IκBs, and the release of active NF-κB that translocates into the nucleus to activate the expression of several genes, including those encoding for the IκB inhibitors (Hoffmann et al., 2002). Re-expression of IκBs contributes to relocate NF-κB in the cytoplasm, which is an inhibitory feedback loop for the system. A second feedback loop is centred on the protein A20 (Ashall et al., 2009), which upon activation of the system inhibits the IKKs. Thus, these two layers provide a nonredundant regulation of NF-κB response to external stimuli (Werner et al., 2008). Of note, the IκB inhibitors also show variable sensitivity to different stimuli and respond on different timescales to fine-tune the signaling pathway (Paszek et al., 2010; Shih et al., 2009). This system of negative feedback loops (summarized in Figure 1A upper panel) provides both a tight control of the response to external stimuli and flexibility. As a consequence of this complex wiring, upon constant stimulation the nuclear concentration of NF-κB in each cell oscillates with heterogeneous dynamics according to each cell’s susceptibility and to the inherent stochasticity of the system (Nelson et al., 2004; Tay et al., 2010; Zambrano et al., 2014a). Due to such asynchrony at the single cell level, the NF-κB response appears almost non-oscillating at the cell population level (Hoffmann et al., 2002), and it is difficult to correlate NF-κB oscillatory profiles with gene expression outputs.

Figure 1 with 11 supplements see all
Periodic forcing turns damped heterogeneous oscillation into synchronous sustained oscillations.

(A) The activity of NF-κB is regulated through different negative feedbacks provided by the inhibitors IκB and A20. The scheme at the bottom represents a generic forcing with periodically alternating TNF-α doses D1 and D2 of duration T2+T1 = Tf; Tf is the period of the forcing. (B, C) Oscillations observed in three GFP-p65 cells obtained by computing the nuclear to the cytoplasmic GFP intensity (NCI) for constant flow of 10 ng/ml TNF-α (B) and upon alternating doses D1=10 ng/ml TNF-α, D2=0 ng/ml and T1=T2 =45 min (C). Each colour corresponds to a single cell trace. Oscillatory patterns can be effectively visualised using the phase ϕ of the oscillation, which is 2π in the maxima of the oscillatory peaks (yellow) and π in the local minima (green) in the colour phase-plot of ϕ(t) below the panels. Scale-bar for ϕ is on the right. (D) Time lapse images of cells under constant stimulation displaying the characteristic heterogeneous nuclear-to-cytoplasmic translocations. (E) Phase plot drawn for 50 cells, of 105 analysed, showing the asynchrony of the oscillations except for the first peak. (F) Distribution of the experimentally computed period of the oscillations Texp, measured as the time between two consecutive oscillatory peaks. The distribution has a maximum at T0 =90 min, which corresponds to the natural period. (G) Quantification of the height for each peak. (H) Time lapse images of the cells under periodic stimulation, showing synchronous NF-κB translocations between cytoplasm and nucleus. (I) Phase plot for 50 cells, of 206 analysed, showing a clear synchrony of the oscillations. (J) Distribution of the period of the oscillations Texp. Texp corresponds almost perfectly to the period of the forcing. (K) Quantification of peaks height variation as described in G; values for n>1 are slightly higher than those observed under constant stimulation. Figure supplements from 1 to 10 are provided.

https://doi.org/10.7554/eLife.09100.003

Obtaining synchronous NF-κB dynamics in cell populations is therefore important to study gene expression, and can give insights into the dynamics at a collective level. White’s group showed that short trains of TNF-α pulses produce rounds of synchronous NF-κB translocation; the stimulation frequency affected both the translocation amplitude and gene expression levels (Ashall et al., 2009). A study performed while the present one was in progress indicated that entrainment of NF-κB oscillations at population level arises when 3T3 cells stably transfected with GFP-p65 are perturbed with regular trains of sawtooth-like profiles of TNF-α (Kellogg and Tay, 2015). Here, we used GFP-p65 mouse embryonic fibroblasts (MEFs) derived from knock-in mice and thus expressing physiological levels of p65 (De Lorenzi et al., 2009). Using a microfluidic device, we exposed cells to well-defined, periodic TNF-α stimuli, and eliminated continuously both catabolites produced by the cells and secreted proteins, which might generate a secondary autocrine/paracrine response. We find that GFP-p65 MEFs lock their NF-κB oscillations to the periodic external signal, and become synchronized. However, the 1:1 locking is maintained over a wide variety of frequencies of the driving TNF-α stimulus, does not improve over repeated stimulation and actually disappears fast when the external stimulus ceases, all of which suggest that entrainment is not achieved. The mathematical model we developed indeed suggests that NF-κB can behave as a damped oscillator, analogous to a mechanical damped harmonic oscillator; damped oscillators do not entrain but follow external forcing while it is present. Taking advantage of the fact that GFP-p65 MEFs under periodic stimulation behave as a synchronous population, we analyzed the transcriptional output at genome-wide level. We find that one single NF-κB dynamics translates into 3 different dynamics of transcriptional regulation, all of which can be reproduced by our mathematical model; the key discriminator is the parameter representing mRNA degradation. Furthermore, the three dynamical patterns correspond to specific functions, which suggests that group of genes were positively selected by evolution.

Results

Periodic forcing turns heterogeneous NF-κB oscillations into synchronous oscillations

To characterize the response of the NF-κB oscillatory system to different external stimuli, we made use of GFP-p65 knock-in MEFs (De Lorenzi et al., 2009; Sung et al., 2009; Zambrano et al., 2014a) cultured in a microfluidic device to control precisely the concentration and timing of the TNF-α stimulation (Materials and methods and Figure 1—figure supplement 1A and B). Notably, the flow rate in the microfluidic device is constant, so that the concentration of TNF-α cannot change due to the activity of the cells, nor can the medium accumulate catabolites or secreted proteins.

Constant stimulation of GFP-p65 knock-in cells in the microfluidic plate with a constant flow of 10 ng/ml TNF-α for up to 15 hr induced nuclear-to-cytoplasmic p65 oscillations (Figure 1B; Video 1 and Figure 1—figure supplement 2D) that are qualitatively similar to the heterogeneous oscillations observed with static stimulation (Sung et al., 2009; Zambrano et al., 2014a).

Several controls were performed to exclude possible damaging effects related to imaging. Cells under constant flow of fresh medium without TNF-α divide and show almost no cell death (Video 2); cell death observed under continuous flow depends on the dose of TNF-α (Figure 1—figure supplement 8 and Video 1) independently of the imaging conditions. UV-induced DNA damage due to imaging (Cadet et al., 2005) is negligible as assessed by immunostaining for thymine dimers (Komatsu et al., 1997; Sinha and Hader, 2002). DNA Damage Response (DDR) is also negligible as assessed by immunostaining of imaged cells for gammaH2AX (Marti et al., 2006; Oh et al., 2011; Staszewski et al., 2008) (Figure 1—figure supplement 5 and 6). Moreover, neither Hoechst exposure nor TNF-α affect dramatically the cell cycle (Figure 1—figure supplement 11). The controls to exclude possible effects of the nuclear dye (Ge et al., 2013; Martin et al., 2005) or photo-damage (Cole, 2014) are described in Materials and methods and in figure captions.

Video 1
Dynamics for constant flow of TNF-α.

Imaging of GFP-p65 knock-in cells in the microfluidic chamber stimulated with a constant flow of 10 ng/ml TNF-α for 12 hr. The stimulus induced nuclear-to-cytoplasmic p65 oscillations that are qualitatively similar to the heterogeneous oscillations observed with static stimulation. Along time, it is possible to appreciate the increase of TNF-induced cell death and apoptosis.

https://doi.org/10.7554/eLife.09100.015
Video 2
Dynamics for Hoechst-stained cells in the absence of TNF-α.

Several controls were performed to exclude possible damaging effects related to imaging. This video shows that cells under constant flow of fresh medium without TNF-α divide and very few events of cell death are detectable. Left part: GFP channel; right part HOE channel. Twelve-hour imaging.

https://doi.org/10.7554/eLife.09100.016

To analyse the oscillatory dynamics, we refined our recently published pipeline (Zambrano et al., 2014a) to compute the nuclear to cytoplasmic intensity (NCI) of GFP-p65 fluorescence for hundreds of cells per condition (Materials and methods and Figure 1—figure supplement 1C). Although this measure implies a partial segmentation of the cytoplasm and thus is more elaborate than the background-corrected mean nuclear intensity used recently by different authors (Kellogg and Tay, 2015; Lee et al., 2014; Lee et al., 2009; Sung et al., 2014) it is a ratio of intensities and thus it self-corrects for changes in image intensity in our experimental settings (see e.g. Video 1 or Video 2). Furthermore, provided that the total amount of p65 is constant (Figure 1—figure supplement 10), NCI(t) is a monotonic function of nuclear NF-κB (Figure 1—figure supplement 9B) so that oscillations are observed in the former if and only if they are present in the latter (further details in Materials and methods). Finally, a rigorous description of the observed dynamics was achieved through the automatic identification of statistically significant peaks (Materials and methods and Figure 1—figure supplement 2A and B).

The timing between two consecutive peaks is the experimental period Texp. Texp has a distribution with a maximum at 90 min (Figure 1F), which is the approximate intrinsic period reported for NF-κB oscillations that we called T0 (Nelson et al., 2004; Tay et al., 2010; Zambrano et al., 2014a). As observed for static stimulation (Zambrano et al., 2014a), the observed oscillations are very heterogeneous (Figure 1B) and asynchronous (Figure 1—figure supplement 3A and Figure 1—figure supplement 2D). However, the variety of behaviors observed, including non-oscillating cells, is compatible with the dynamics observed in the same cells for static culture conditions, see (Sung et al., 2009; Zambrano et al., 2014a). Similar behaviors are found in the absence of nuclear staining (see Figure 1—figure supplement 9C), which led us to conclude that imaging conditions do not interfere with NF-κB signaling.

The tail in the period distribution (Figure 1F) and the observed average of four peaks in 12 hr of stimulation (Figure 1—figure supplement 3B), instead of the expected eight, indicate that in our cells the heterogeneous oscillations are damped and tend to converge to an equilibrium state under stimulation (Figure 1B). Indeed, although oscillatory peaks are observed for most of the cells, they are infrequent and irregular for times beyond 6 hr (Figure 1B and Figure 1—figure supplement 2D) as previously reported for the same cells in static conditions (Sung et al., 2009; Zambrano et al., 2014a). Hence, we describe here these heterogeneous oscillations as damped in contrast with the sustained oscillations, which continue regularly and unabated for a very long time in continuously stimulated cells, see for example (Kellogg and Tay, 2015).

Damped oscillations can easily emerge in the NF-κB genetic circuit. Indeed, a minimal deterministic model that takes into account the basic elements of the NF-κB genetic circuit (Zambrano et al., 2014b) (Figure 1A, Figure 1—figure supplement 4, see Materials and methods for a complete description of the model) shows that different combinations of the parameters can lead to different dynamics. The model parameters that we used for our explorations are provided in Supplementary file 2. We denote as PS those specifying the external signal and as PNF-κB those used to model the double IκB and A20 negative feedback; in our explorations, we allow them to vary differently depending on the associated uncertainty about their values (Materials and methods and Supplementary file 2). We generated a library of randomized parameters and found that the system presents a fixed point whose stability changes depending on the parameters (details on the stability analysis are found in the Materials and methods section). The vast majority of parameter combinations give rise to damped oscillations (when all the eigenvalues have all negative real parts, see Figure 1—figure supplement 5A,B). A smaller fraction of parameter combinations give rise to trajectories that converge to a stable limit cycle around the unstable fixed points (so certain eigenvalues have positive real parts, see Figure 1—figure supplement 5A,B). Interestingly, parameters for oscillating and non-oscillating cells are in similar intervals suggesting that it is the precise combination of the parameters, rather than a single one, what determines the resulting dynamics (see Figure 1—figure supplement 5C). Our simulations might also explain why other researchers found continuous periodic oscillations with T0 = 90 min under a constant flow of TNF-α (Kellogg and Tay, 2015, and see Discussion) and the variety of damped oscillatory dynamics upon LPS recently reported for fibroblasts (Cheng et al., 2015) and macrophages (Sung et al., 2014). Our exploration shows further how variations of the parameters can give rise to a variety of dynamics that reflects what we find in an isogenic population (Figure 1B and Figure 1—figure supplement 2).

Considering the heterogeneity of dynamics, to better visualize the collective oscillatory state of the population in each condition, following Mondragon-Palomino et al., 2011, we computed and represented the phase of the oscillation ϕ(t) for each cell by detecting peaks and setting ϕ=0 (2π) at the maximum of each peak and ϕ(t)=π in the minimum between two peaks (phase plots for the green time series are depicted in Figures 1B,C. See also Materials and methods and Figure 1—figure supplement 2B,C). Time series for single cells (Figure 1B and C) were converted to phase plots (Figure 1E,I) where each row represents one cell. Thus, oscillatory peaks can be easily observed. In the phase plot, the first response to constant TNF-α right after t=0 hr is synchronous in the population, but this synchrony is quickly lost, as previously reported (Nelson et al., 2004; Tay et al., 2010; Zambrano et al., 2014a).

To investigate the response of the NF-κB oscillator to perturbations in the cell’s environment, we switched periodically the stimulus concentration in the culture chambers in the microfluidics apparatus. TNF-α switching between doses D1 and D2 occurs in less than 1 min and generates a tightly controlled square profile of stimulation. Stimuli were applied for intervals of time T1 and T2. We refer to Tf = T1+T2 as the period of the forcing and to D1–D2 as the amplitude of the forcing (Figure 1A, lower panel).

We started our analysis by applying a periodic stimulation of 90 min, which is close to the intrinsic period of the NF-κB oscillatory system. Single-cell traces are provided in Figure 1C for cells stimulated with D1 =10 ng/ml TNF-α, D2 =0 ng/ml with Tf=90 min. (T1=45 min and T2 =45). Peaks are present in all the forcing cycles (Figure 1C and Figure 1—figure supplement 3C) and synchronous, in contrast to the asynchronous peaks observed in constantly stimulated cells (Video 3). Hence, we observe in this condition – a square forcing – a forcing-induced synchronous dynamics reminiscent of that obtained by applying short trains of pulses (Ashall et al., 2009). Our synchronous oscillations are visually apparent both in time-lapse images (Figure 1H) and phase plots (Figure 1I). Texp is sharply distributed around 90 min (Figure 1J), corresponding to the expected eight peaks in 12 hr (Figure 1—figure supplement 3D). The height of the first peak (Figure 1K) is similar to the height of the first peak under constant stimulation (Figure 1G), indicating that the experimental conditions are highly comparable. The height of the subsequent peaks is slightly higher for periodically stimulated cells but still heterogeneous.

Video 3
Dynamics for Tf=90 min.

Imaging of cells stimulated with D1 =10 ng/ml TNF-α, D2 =0 ng/ml and Tf=90 min. (T1=45 min and T2 =45). Peaks are present in all the forcing cycles and synchronous. Twelve-hour imaging.

https://doi.org/10.7554/eLife.09100.017

Taken together, these results indicate that a periodic external stimulus can lock in step a population of cells in a long series of synchronous NF-κB oscillations, in clear contrast with the remarkable oscillatory asynchrony and heterogeneity under constant stimulation.

Oscillatory synchrony increases with the forcing amplitude

We then systematically analysed the response to variable forcing amplitudes, by applying different concentrations of TNF-α. The phase plots obtained for TNF-α ranging from 10 to 0.1 ng/ml with D2 = 0 ng/ml show that the coherence of the oscillations decreases with decreasing doses of TNF-α but still persists at low doses (Figure 2A and Figure 2—figure supplement 1). Texp is less sharply distributed around 90 min as the forcing amplitude decreases. Concomitantly, a second peak corresponding to multiples of Tf becomes progressively more conspicuous (Figure 2B), suggesting that for lower doses of TNF-α an increasingly larger fraction of cells do not respond with a peak to some of the forcing cycles. In the plots of peak maxima (Figure 2C) the height of the first peak is not reproduced by the second and the third peaks, compatible with what reported for this forcing period (Ashall et al., 2009). However the heights converge under repeated forcing to a constant value for all the forcing amplitudes. This is typical of nonlinear dissipative oscillators in the presence of periodic forcing (Goldstein et al., 2001).

Figure 2 with 3 supplements see all
Synchronous oscillations arise for different forcing amplitudes.

(A) Representative phase plots for 25 cells (out of 216, 80, 188, 263, 225 cells analysed) stimulated with D1=10, 1, 0.5, 0.25, 0.1 ng/ml TNF-α, D2=0 ng/ml, for Tf =90 min, T1=T2=45 min. (B) Distributions of the periods for the cells shown in panel (A); distributions become narrower as the dose D1 increases. The appearance of a second peak at Texp=3 hr at lower doses means that in some cycles a fraction of cells miss a peak and the interval to the next one is double. (C) Quantification of height of the nth peak in the different conditions considered in (A). By decreasing the stimulus amplitude, the ratio tends to stabilize to a constant value. (D) Distribution of the phase difference Δϕ for the forcings considered: Δϕ becomes narrower as the forcing amplitude is increased. (E) The synchrony intensity η, an entropy-based measure on how widely distributed the values of Δϕ are, increases with the amplitude of D1 for Tf =90 min. (F) The synchrony intensity ηn computed using only the peaks observed in each forcing cycle shows that the synchrony does not increase with the successive cycles of forcing. Figure supplements from 1 to 3 are provided.

https://doi.org/10.7554/eLife.09100.018

The degree of synchrony for each experimental setting of forcing can be evaluated considering the distribution of the phase difference Δϕ between the timing of each oscillatory peak and the beginning of the forcing (see Materials and methods). For asynchronously oscillating cells under constant TNF-α stimulation, such distribution is flat (Figure 2—figure supplement 2). When cells are stimulated with different forcing amplitudes, the distribution of Δϕ is narrower for higher doses of TNF-α (Figure 2D), indicating a higher degree of locking of the oscillations to the forcing for increasing amplitudes.

We then compared the degree of synchrony of the cell populations using the synchrony intensity η, an entropy-based quantifier of the distribution of the phase difference Δϕ that is 0 for flat distributions and 1 for delta-like distributions (Mondragon-Palomino et al., 2011, and Materials and methods). Synchrony intensity increases with the dose, but is nonzero even when D1=0.1 ng/ml TNF-α (Figure 2E). To understand if D2=0 was a necessary condition for synchrony we applied D2>0. Results show that the synchrony is maintained when D1 is threefold D2, while for smaller differences synchrony is almost lost (Figure 2—figure supplement 3).

We next investigated whether synchrony improved under repeated forcing. We computed ηn, the synchrony intensity at the n-th cycle of forcing (i.e. Δϕ computed for peaks in the time intervals [(n–1)Tf,nTf) for n≥1). We find that ηn does not increase (Figure 2F) in successive cycles of external forcing, meaning that the oscillations do not become more synchronous as the system is perturbed repeatedly.

Taken together, the above results indicate that NF-κB dynamics adapts to varying amplitudes of periodic inputs. However the system does not seem to learn from the previous periodic forcing cycles.

Cells oscillate synchronously following a variety of forcing amplitudes and periods

We then tested the adaptability of the NF-κB system to different forcing periods Tf, ranging from 0.5T0 to 2T0.

We first considered periodic perturbation of Tf=2T0=180 min, with T1=30 min and T2 =150 min, and TNF-α doses D1 ranging from 10 ng/ml to 0.1 ng/ml, D2 =0 ng/ml. Oscillations are locked to the forcing even for D1 =0.1 ng/ml (Figure 3A, see also Figure 3—figure supplement 1,A–C and Video 4). The use of T1=30 min assures the existence of sharp oscillatory and transcription peaks (see below) and also leads to synchrony of the oscillations for Tf=T0 (see Figure 3—figure supplement 3, bottom panels). For all doses we find a bimodal distribution of the oscillation periods, with an overall maximum at 180 min and a much smaller relative maximum at 90 min (Figure 3C). This indicates that after a single pulse of forcing some cells oscillate a second time with a period similar to the intrinsic one. However, a closer analysis of peak maxima for time intervals of the form [(n–1)T0, nT0) = [(n–1)Tf/2, nTf/2) (Figure 3E) reveals that peaks arising right after each periodic stimulation (even n) are higher than the rare ones (odd n) arising when stimulation is absent.

Figure 3 with 3 supplements see all
Cells adjust oscillations to different periods for a wide range of forcing amplitudes.

Representative phase plots for 25 cells stimulated with D1=10, 1, 0.1 ng/ml TNF-α (out of 151, 77, 123 cells analysed, respectively) D2=0 ng/ml for (A) Tf =180 min with T1= 30 min and (B) Tf =45 min with T1=22.5 min (analysed 101, 112, 119 cells). (C, D) Plots showing the distributions of the periods for the conditions given in panels (A) and (B), respectively. (E, F) Average peaks height in the intervals [(n–1) T0, T0) and [(n–1) T0/2, nT0/2), for A and B, respectively. In (E), even n correspond to peaks right after stimulation, odd n correspond to the small peak arising between two consecutive stimulations. (G, H) Distribution of the phase difference Δϕ for the forcings in A and B: Δϕ has narrower distributions for higher doses. (I) The synchrony intensity η grows with the doses for Tf=180 min (blue) and Tf=45 min (red). (J, K) Synchrony intensity plots show that ηn does not increase as successive cycles of forcing are applied to the system, both for Tf=180 min and Tf=45 min, respectively. All the analyses included all the tracked cells with no preselection of the responding ones. Figure supplements 1 to 3 are provided.

https://doi.org/10.7554/eLife.09100.022
Video 4
Dynamics for Tf=180 min.

Imaging of cells stimulated with D1 =10 ng/ml TNF-α, D2 =0 ng/ml and Tf=2T0=180 min (T1=30 min and T2 =150). Oscillations are locked to the forcing. Twelve-hour imaging.

https://doi.org/10.7554/eLife.09100.026

Oscillations in each condition show a good degree of synchrony, as described by the distributions of the phase differences Δϕ, which are remarkably narrow for all doses (Figure 3G); the synchrony intensity η is nonzero even for D1 = 0.1 ng/ml and increases with the dose (Figure 3I, blue line). When D2>0 ng/ml we found synchronized oscillations for D1 at least three times D2 (Figure 3—figure supplement 2).

We considered also the ability of the NF-κB system to respond to stimulations of periodicity below the intrinsic one. We selected Tf=T0/2=45 min, but T1=30 min led to a poor synchrony (see Figure 3—figure supplement 3, top panels). We therefore reduced T1 to 22.5 min, and thus T2=22.5, which proved to be a sufficiently long resetting of the external signal; in these conditions we obtained a sharply defined dynamical response. Oscillations are locked in step with D1=10 and 1 ng/ml (Figure 3B, Video 5), whereas synchronization is weaker for D1=0.1 ng/ml (this is also evident in the average NCI dynamics, see Figure 3—figure supplement 1,D–F). The experimentally determined oscillatory period Texp (Figure 3D) is narrowly distributed around Tf=45 min for D1=10 ng/ml, but we also find a small peak around Texp=90 min, indicating that some cells skip the oscillation elicited by the forcing. The peak height for each forcing cycle tends to stabilize to a constant value and to be lower for lower forcing amplitude, although the differences are small (Figure 3F). Also in this case, the synchrony correlates with the dose: in fact, the phase difference Δϕ tends to be more narrowly distributed for higher values of D1 (Figure 3H), and the synchrony intensity η grows with the dose (Figure 3I, red line).

Video 5
Dynamics for Tf=45 min.

Imaging of cells stimulated with D1 =10 ng/ml TNF-α, D2 =0 ng/ml and Tf=T0/2=45 min (T1=22.5 min and T2 =22.5 min). Such short wash-out provides a sufficient resetting of the external signal; in these conditions we obtained a sharply defined dynamical response and oscillations are locked in step. Twelve-hour imaging.

https://doi.org/10.7554/eLife.09100.027

It has been recently reported that NF-κB oscillations can be synchronized by entrainment (Kellogg and Tay, 2015). However, in our experiments the synchrony intensity for each cycle, ηn, does not increase when the system is forced repeatedly, either with periods of 90, 180 or 45 min (Figures 2F, 3J,K), in contrast to what is observed for entrained oscillators (Mondragon-Palomino et al., 2011). Furthermore, when oscillators are entrained by external forcing, their oscillations also show m:n resonant patterns (m oscillations for each n cycles of the forcing) in so-called Arnold tongues (Pikovsky et al., 2003). However, we only observed 1:1 synchronization patterns for all the periods of forcing (Figure 3A,B), with no clear evidence of the 2:1 and 1:2 patterns expected for entrainment. We thus investigated further whether the synchronization mechanism we observed does correspond to entrainment.

The synchronization mechanism of GFP-p65 knock-in MEFs is not entrainment

Once the forcing ceases, entrained oscillators dephase gradually, losing their sharply defined common entrained oscillatory period. We then investigated whether our cells kept a memory of the synchronous dynamics once the periodic forcing ceased. The phase plots for 75 cells forced with Tf=45 min D1 =10 ng/ml (Figure 4A, washing out after 16 forcing cycles), and for Tf=90 min with D1 =1 ng/ml (Figure 4B, washing out after eight forcing cycles), indicate that when the external stimulus stops, cells lose quickly their synchrony. This is also evident from single-cell NCI traces and from the peaks detected in the same conditions (Figure 4C and D respectively). Indeed, only for Tf=45 min we observed small peaks after the last cycle of the forcing, but these peaks are on average 90 min away from the last forced peak, and correspond to the intrinsic oscillatory period of our oscillator. The quick loss of synchrony can also be observed in the evolution of the synchrony quantifier ηn for the two last cycles of the forcing (Figure 4E, n=1,2) and the two subsequent (virtual) ones (Figure 4E, n=3,4). This is also the case when considering cells forced with Tf=90 min and a high stimulus amplitude D1 =10 ng/ml, after which we apply a constant flow of 10 ng/ml TNF-α: cells present a last well-defined translocation (Figure 4F and 4G) and rapidly dephase (Figure 4H).

Figure 4 with 4 supplements see all
Cells do not keep a memory of the synchronous oscillatory dynamics.

(A, B) Phase plots of the last two oscillation cycles for Tf =45 min (D1=10 ng/ml) and Tf =90 min forcing (D1= 1 ng/ml) (number of forcing cycles are indicated above) followed by a period of 3 hr with no stimulation (75 cells are displayed out of 106 and 197 cells analysed, respectively). (C, D) NCI time series at single cell level (green lines) for the two conditions (A) and (B). Blue triangles indicate the peaks considered in the computing. The thick black line is the average NCI, showing a small peak 90 min after the last forced peak for Tf =45 min (C), compatible with the natural timescale of the free oscillations. (E) The synchrony intensity ηn for the last two forced peaks (n=1, 2) and for the peaks detected in the absence of the stimulus (n=3, 4) illustrates the fast loss of synchrony. (F) Phase plots of the last two oscillation cycles for Tf =90 min (number of forcing cycles are indicated above) (D1=10 ng/ml) followed by 4.5 hr flow of 10 ng/ml TNF-α (75 cells are displayed out of 149 cells analysed). (G) NCI time series at single-cell level (green lines). Blue triangles indicate the peaks considered in the computing. The thick black line is the average NCI, showing a small peak 90 min after the last forced peak, compatible with the natural timescale of the free oscillations. (H) The synchrony intensity ηn for the last two forced peaks (n=1, 2) and for the peaks detected in the presence of the stimulus (n=3,4) illustrates the fast loss of synchrony. Figure supplements from 1 to 4 are provided.

https://doi.org/10.7554/eLife.09100.028

Kellogg and Tay (2015) found no synchrony for Tf=60 min. However, we found that our cells can actually be synchronized under periodic stimulations with Tf=60 min (T1 =T2 =30 min) when D1 =10 ng/ml (Video 6) and D2 =0 ng/ml (Figure 4—figure supplement 1A). Furthermore, we found that synchrony is still visible for D2>0 ng/ml (Figure 4—figure supplement 1B–D), as we had observed for other stimulation periods. The difference might arise from the fact that we use a forcing different from the sawtooth-like TNF-α profile used by Kellogg and Tay (2015), which is obtained by periodically and quickly replacing the medium in contact with the cells with fresh TNF-α. The sawtooth concentration profile is assumed to arise due to a clearance process of TNF-α by degradation and internalization. We then applied sawtooth forcing to our cells (Figure 4—figure supplement 4). Cells lock their oscillations to the sawtooth forcing of periods Tf=60 min and Tf=90 min at the doses considered. Interestingly, the synchronization is lost for Tf=180, suggesting that autocrine-paracrine signalling might introduce a distortion in the cell environment, which is reduced when the medium is changed more frequently. Alternatively, the different geometry of our microfluidic with respect to the one used in Kellogg and Tay (2015) might lead to a slower TNF-α degradation and produce profiles closer to static concentration than to oscillatory stimulation. The profiles observed are indeed similar to some of the dynamics observed for alternating doses of TNF-α with D2>0 ng/ml, such as those for Tf=60 min (Figure 4—figure supplement 1B–D), Tf=90 min (Figure 2—figure supplement 3) and Tf=180 min (Figure 3—figure supplement 2), with compatible values of η. However, it is clear that with sawtooth forcing of Tf=60 min the synchrony is strong, at variance with the results reported by Kellogg and Tay (2015).

Video 6
Dynamics for Tf=60 min.

GFP-p65 cells can be synchronized under periodic stimulations with Tf=60 min (T1 =T2 =30 min) when D1 =10 ng/ml and D2 =0 ng/ml. Twelve-hour imaging.

https://doi.org/10.7554/eLife.09100.033

Overall, we conclude that the synchronization mechanism that we observe is not entrainment. Rather, it is similar to that of simple mechanical damped oscillators, such as the damped harmonic oscillator, which after a perturbation tend to relax to an equilibrium state. When these mechanical systems are challenged with an external periodic forcing, the period of the oscillations tends to match the period of the forcing (Pikovsky et al., 2003). In fact, our minimal mathematical model of the NF-κB circuit (Figure 1—figure supplement 4) reproduces damped oscillations (e.g. Figure 1—figure supplement 5A,B,D). In these conditions, the period of the periodically forced system converges to the period of the forcing (Figure 4—figure supplement 2). Of note, the model reproduces the small peaks appearing between two forcing cycles (Figure 4—figure supplement 3), similar to the small peaks that we observed for the same forcing as in Figure 3A.

When behaving as a damped oscillator, the NF-κB system is well suited to quickly synchronize its oscillatory period to input signals of a wide variety of timescales. This is probably a desired feature for a system that underlies responses to environmental challenges, which do not come with a particular periodicity. This is the reverse of what might be desirable for circadian clocks, whose design should privilege the entrainment to the 24-hr light/darkness period. To understand how the oscillatory pattern of NF-κB affects its biological output, we focused on the transcription dynamics of genes controlled by NF-κB.

Synchronous NF-κB dynamics reveal different patterns of gene expression

The heterogeneity of NF-κB dynamics in the cell population under constant stimulation makes it difficult to establish a direct correspondence between NF-κB activity and transcription. However, under conditions of synchronous oscillation the average dynamics corresponds to the dynamics of single cells. This is for example the case for cells stimulated with D1 =10 ng/ml D2=0 ng/ml TNF-α, T1=30 min and T2 =150 min (Tf= 2T0) (Figure 5A). Numerical simulations with a simple mathematical model of transcription under the control of the regulatory network using parameters from our previous work (Figure 5—figure supplement 1, Materials and methods and Zambrano et al., 2014b) suggested that for some genes there should be a clear coordination between the input pulsed signal, the oscillating dynamics of NF-κB and the transcriptional output (Figure 5B).

Figure 5 with 2 supplements see all
Synchronous NF-κB oscillations lead to population-level coordinated transcription.

(A) NCI plot of single cell oscillations (green lines) and population average (black lines) for cells stimulated with D1=10 ng/ml TNF-α, D2=0 ng/ml, T1=30 min and T2=150 min. The open circles represent the fitting obtained using our minimal mathematical model. (B) The mathematical model predicts waves of transcription (orange plot, right) coordinated with the stimulus (red plot, top) and p65 oscillations (green plot, left). (C, D) q-PCR time course of nascent and mature mRNAs for the prototypical early and late genes IκBα and Ccl5, respectively. (E, F) Transcription profiles for mature IκBα (red) and Ccl5 (blue) RNAs (dots) can be accurately fitted (lines) by our minimal mechanistic mathematical model. The fittings were performed by keeping common the parameters regulating the external signal (PS) and the dynamics (PNF-κB) in (A), (E) and (F) but using different gene expression parameters PG for (E) and (F). Figure supplements 1 to 2 are provided.

https://doi.org/10.7554/eLife.09100.034

To validate the model’s prediction, we tested RNA transcription in the synchronous oscillatory condition shown in Figure 5A. RNA samples were prepared at time 0 and 20, 40 and 60 min after each TNF-α pulse. Quantitative RT-PCR was performed for both the nascent unspliced and the mature mRNA forms of the prototypical early and late genes IkBa and Ccl5/RANTES, respectively (Figure 5—figure supplement 2 and Materials and methods). Nascent transcription of both genes starts synchronously with each TNF-α pulse and p65 nuclear localization (Figure 5C,D), and the levels of nascent transcripts clearly oscillate. Mature transcripts of the early gene IκBα follow the oscillatory dynamics, show transcription-coordinated splicing and do not accumulate after stimulus wash-out (Figure 5C). On the contrary, the mature form of the late gene Ccl5 accumulates slowly and progressively along 12 hr (Figure 5D).

This latter observation extends the notion that the splicing rate for inflammatory genes plays a key role in regulating the timing of gene expression, as previously reported (Hao and Baltimore, 2013). The induction levels for IκBα in the first hours hardly differ from those under chronic stimulation ((Hao and Baltimore, 2013); Figure 5 and Figure 5—figure supplement 2). Interestingly, we also find oscillatory dynamics of IκBα and Ccl5 nascent transcripts when cells are stimulated with a forcing period of Tf= 90 min, although the oscillations are less prominent than for Tf=180 min. The expression of Ccl5 for Tf= 90 min grows monotonically (Figure 7—figure supplement 1).

We then asked whether our minimal mathematical model (Figure 5—figure supplement 1 and Figure 1—figure supplement 4) could simultaneously fit NCI dynamics (Figure 5A) and transcription profiles. Our model of transcription can be interpreted as a population-level NF-κB–driven telegraph model (Suter et al., 2011) that includes non-cooperative transcription activation by NF-κB (Giorgetti et al., 2010) and the role as transcriptional repressor reported for IκBα (Arenzana-Seisdedos et al., 1995; Kellogg and Tay, 2015; Lipniacki et al., 2004; Lipniacki et al., 2006a; Lipniacki et al., 2006b; Tay et al., 2010); we denote the parameters used to fit the transcription as P(details in the Materials and methods section). The fittings were performed by keeping constant the parameters of the external signal PS (which is externally imposed) and using the same parameters regulating NF-κB dynamics, PNF-κB in Figure 5A,E,F, while using different individual transcription parameters PG for each gene (Figure 5E,F). Despite being much simpler than other existing models of NF-κB dynamics (Ashall et al., 2009; Tay et al., 2010), our model fits faithfully the observed transcription dynamics for both IκBα (Figure 5E) and Ccl5 (Figure 5F).

Overall, our results show that the expression of two different genes can follow different dynamical patterns even if the entire cell population is effectively locked to the dynamics of the transcription factor that controls both genes.

Transcription dynamics discriminates between groups of functionally related genes

Genome-wide profiling represents the logical further step to understand whether additional NF-κB regulated genes follow the diversified dynamics observed for Nfkbia/IκBα and Ccl5. We performed microarray analysis on RNA purified from GFP-p65 knock-in cells under either constant or pulsed TNF-α stimulation (Zambrano et al., 2016). Transcriptional outputs were subjected to unsupervised clustering of standardized profiles to group the transcription profiles by their shape, while minimizing differences in fold-changes (Materials and methods). For constantly stimulated cells, gene expression profiles reproduce well the previously observed kinetics (Rabani et al., 2011; Sivriver et al., 2011), in which the maximum expression level is achieved with different timings for different genes, and no evidence of the oscillatory dynamics of NF-κB is discerned (Figure 6—figure supplement 1).

We then quantified genome-wide expression profiles in the population shown in Figure 5, which oscillates synchronously in the forcing regime of Tf=180 min (Materials and methods). We found 970 genes whose expression responds to TNF-α stimulation, of which 499 increase and 471 decrease relative to unstimulated cells. Genes distributed in six highly homogeneous clusters; a higher number of clusters did not reduce significantly the inter-cluster distance (Figure 6—figure supplement 2A). Three of the clusters contain genes with increasing expression and three contain genes with decreasing expression (Figure 6A, Figure 6—figure supplement 2). The clusters display profiles that range from the oscillatory dynamics of IκBα (Cluster 1) to the slowly increasing dynamics of Ccl5 (Cluster 3). Notably, while genes contained in the three clusters of increasing expression are enriched with known NF-κB target genes (10-14 < p < 0.05, Fisher’s exact test, Materials and methods and Figure 6—figure supplement 3), clusters with decreasing expression are not significantly enriched. The same is observed in chronic stimulation (Figure 6—figure supplement 4).

Figure 6 with 7 supplements see all
Synchronized NF-κB dynamics translates into functionally different dynamical patterns of gene expression, each corresponding to distinct pathways.

(A) Clusters 1–6 were obtained by an unsupervised k-means-like clustering from the genome-wide transcription profiling of samples harvested in the experiment shown in Figure 5A. Line colours are indicative of the membership value of each gene (colour scale at the bottom). Three clusters contain genes with increasing expression (1–3) and three with decreasing expression (4–6). On the y-axis, standardized expression profile in arbitrary units (see Materials and methods). (B) Plots show single-gene mRNA traces (median: thick blue line; 85% and 15% intervals, thin blue lines). The time courses can also be fitted using our minimal mathematical model: shown is the median of the single-gene fits (thick red line) and the 85% and 15% intervals (thin red lines). Fittings were performed using the same parameters for the external signal (PS) and the dynamics (PNF-κB) as in Figure 5, but using different gene expression parameters PG for each gene. (C) Top five pathways of hierarchical level 2 and 3 in the Reactome database significantly enriched in each dynamical cluster. Dot sizes are proportional to the percentage of genes in the cluster belonging to that pathway. Dot colours identify the corresponding p-values (p-value < 0.05 is set as threshold). Scale bars on the right. (D, E) Heatmaps shows the degree of overlap at gene level (D) and pathway level (E) between each of the 6 clusters. Colour scale bar on the right. Figure supplements from 1 to 6 are provided.

https://doi.org/10.7554/eLife.09100.037

We tested if our model would also reproduce the patterns of gene expression dynamics at a genome-wide level in these conditions, keeping the same parameters PS and PNF-κB as in Figure 5 but using different gene expression parameters PG for each gene. Indeed, the model fits well the dynamical patterns observed experimentally in clusters with increasing expression (Figure 6B, grey lines). The median and the intervals of observed and simulated single gene expression levels match remarkably well (Figure 6B, red and blue lines respectively, individual fittings are shown in Figure 6—figure supplement 6, bottom-left panels); the average relative error per timepoint is less than 10% (Figure 6—figure supplement 5, left, further details on the fittings and metrics used are given in Materials and methods).

Genes in the clusters with decreasing expression cannot be fitted with the same accuracy (the average relative error per timepoint is above 10%, see Figure 6—figure supplement 5, right, individual gene fitting examples shown in Figure 6—figure supplement 6, bottom-right panels), further suggesting that these genes are not controlled directly by NF-κB. Indeed, an analysis of the parameters show that gene expression is essentially fitted for those genes as a simple RNA degradation with nearly no contribution of the gene activation/inactivation processes (low Kon,G and Koff,G, see Figure 5—figure supplement 1 and Materials and methods).

Interestingly, our mathematical model indicates that for a given stimulus impinging on NF-κB, the degradation rate of the mRNA is the key parameter for producing the most distinct dynamical profiles. Of note the highest degradation rates are specific for oscillating genes, while those for slowly increasing genes are two orders of magnitude lower (Figure 6—figure supplement 6, top panels). This further underlines the importance of mature RNA processing and degradation in the definition of different patterns of gene expression.

Although the role of IκBα as a repressor has been proposed in a number of theoretical analyses (Kellogg and Tay, 2015; Lipniacki et al., 2004; Lipniacki et al., 2006a; Lipniacki et al., 2006b; Tay et al., 2010) only limited evidence of this role is experimentally available (Arenzana-Seisdedos et al., 1995). Thus, to investigate the role of IκBα in patterning gene expression, we built a model where IκBα does not act as a direct repressor (Figure 6—figure supplement 7A,B). This alternative model fits the gene expression dynamics of each cluster well (Figure 6—figure supplement 7E), and again better for genes with increasing transcription than for those decreasing (Figure 6—figure supplement 7C). The mRNA degradation rate is again the key parameter to discriminate between the different patterns of gene expression (Figure 6—figure supplement 7D). While our fittings do not provide conclusive evidence on the role of IκBα as a transcriptional repressor, they suggest that this is not needed.

Finally, to understand if the different dynamics observed with Tf=90 min were also present in other conditions, we used the same clustering approach for gene expression profiles of cells stimulated with Tf=90 min, D1 =10 ng/ml and D2 = 0 ng/ml (as in Figure 1C). Oscillations in transcript levels were observed, although they were not as conspicuous as for Tf=180 min (Figure 7A), probably due to the fact that 90 min is not long enough for most transcripts to degrade. Enrichment of target genes is similar to the two previous cases (Figure 7C,D).

Figure 7 with 1 supplement see all
Genome-wide clustering for Tf=90 min.

(A) Clusters obtained by unsupervised k-means clustering from the genome-wide transcription profiling of cells perturbed with a forcing of 90 min, D1=10 and D2=0 ng/ml of TNF-α. Lines colours: blue and red indicate low and high membership values, respectively. (B) Top five pathways of hierarchical level 2 and 3 in the Reactome database significantly enriched in each dynamical cluster. Dot sizes are proportional to the percentage of genes in the cluster belonging to that pathway. Dot colours identify the corresponding p-values (p-value<0.05 is set as threshold). Scale bars on the right. (C, D) Enrichment analysis for NF-κB targets from the clusters displayed in panel A. Two lists of NF-κB targets were considered: left: Gilmore’s web-site (www.bu.edu/nf-kb/); right: data from Brasier and Kudlicki groups (Li et al., 2014). Significance is shown as -Log(p-value), a dashed line marks the threshold of significance at p=0.05. (E, F) Heatmaps show the degree of overlap at gene level (E) and pathway level (F) between each of the 6 clusters. (G) Overlap at a gene level between clusters obtained for 180 min and for 90 min. It is particularly high for Clusters 1, those of oscillating genes. Colour scale bar on the right. Figure supplement 1 is provided.

https://doi.org/10.7554/eLife.09100.045

We then performed a pathway analysis to determine if each of our dynamical clusters were enriched in functionally related pathways (see Materials and methods). Indeed, this is the case (Figure 6C): Cluster 1 is broadly related to chemokines and chemokine receptors and contains so-called early genes, Cluster 2 to the immune system and so-called intermediate genes, and Cluster 3 to extracellular matrix rearrangements and so-called late genes (Rabani et al., 2011; Tian et al., 2005a; 2005b). Thus, the clustering also indicates a precedence order in the articulation of the cell's response. Genes with decreasing expression comprise pathways related to metabolism (Cluster 4), cell cycle (Cluster 5) and chromosome maintenance (Cluster 6).

The overlap between clusters is minimal both at gene and pathway level (Figure 6D and Figure 6E, respectively, and Materials and methods). The correspondence between dynamics and gene function is also fulfilled in cells forced with Tf=90 min (Figure 7B), with a low intercluster overlap at gene level (Figure 7E), and slightly higher overlap at pathway level (Figure 7F). Interestingly, there is a good correspondence between clusters identified upon both 180 and 90 min forcing. In particular Cluster 1 (the oscillating cluster) contains almost exactly the same genes in both conditions (Figure 7G).

Taken together, the above results indicate that different gene expression dynamics identify functionally related categories of genes, suggesting a strict relationship between dynamics and function in the cell's response to external stimuli.

Discussion

The NF-κB system synchronizes to external periodic forcing as a damped oscillator

Our results show that single cells activate NF-κB signalling synchronously to external TNF-α stimulation, and that repeated forcing elicits synchronous NF-κB oscillations at population level. However, NF-κB does not behave as a free oscillator: no cell oscillates constantly in the presence of continuous stimulation (Figure 8A and B left, green lines). Furthermore, synchronization among cells does not improve upon repeated forcing (“training”), and single cells trained for a dozen forcing cycles stop oscillating and dephase (“forget”) as fast as cells stimulated only once. Finally, cells can be synchronized with a wide variety of forcing periods and forcing amplitudes, and synchronize to the external forcing (Figure 8A and B right, green lines) without showing any preference for periods resonating with the NF-κB intrinsic period of 90 min; synchronization is more pronounced with high than with low stimulation amplitude.

NF-κB behaves as a damped oscillator that can synchronize to time-varying external stimuli to produce functionally related transcriptional outputs.

(A) The NF-κB system is able to provide different responses to different inputs, from constant (left) to time-varying ones (right). (B) Our cells show damped oscillations to a constant stimulus (left, green lines), although for other cell types sustained oscillations might be possible (magenta line). Damped oscillations can adapt to timevarying inputs (right) and give rise to synchronous oscillations. (C) These synchronous oscillations produce different patterns of gene expression, from oscillating (left, orange lines) to slowly increasing (right, blue lines) and intermediate dynamics (pink lines, centre). We find that each kind of dynamics is typical for genes involved in different cellular functions.

https://doi.org/10.7554/eLife.09100.047

The ability of NF-κB to oscillate in tune with a forcing of 45 min illustrates well the plasticity of synchronization: 45 min is equivalent to one half of the intrinsic period, and entrained cells would skip one forcing period out of 2 (following the 1:2 resonance), but our data do not support this interpretation, since cells only infrequently skip a forcing period. Simulations with our mathematical model can actually reproduce the ability of the system to synchronize to different forcing periods. The ability of our cells to synchronize to a 60 min periodic forcing points in the same direction.

Thus, we find that our GFP-p65 knock-in cells are not entrained, as do not satisfy the essential precondition – sustained oscillations – and two critical tests for entrainment – increasing synchronization upon repeated forcing, and synchronization to the natural frequency of the internal oscillator (Pikovsky et al., 2003). Our cells appear closer to the classical textbook example of a damped harmonic oscillator (Goldstein et al., 2001), whose frequency corresponds to the frequency of the external stimulation, and whose synchronization to the external stimulus increase monotonically with stimulus amplitude.

We interpret the damped oscillator behaviour of the NF-κB system as an intrinsic characteristic of the NF-κB system which allows it to adapt to a wide variety of inputs and to quickly reset. In fact, pathogens and inflammatory signals do not come in regular periodic patterns, and a carefully calibrated inflammatory response that would adapt itself to any irregular pattern of stimulation would be positively selected for.

We explored our minimal mathematical model to understand why we observed damped oscillations whereas Kellogg and Tay observed sustained oscillations under constant stimulation, a behaviour that presumably is the root of the different synchronization mechanisms observed. Our simulations suggest that sustained oscillations would not be the norm, but rather the result of a specific set of parameters (see Figure 1—figure supplement 5A,B). Indeed, the NF-κB network has an equilibrium that depending on the combinations of parameters can be stable, giving rise to damped oscillations that converge to it, or unstable, giving rise to a stable cycle around it (see e.g. Figure 8B, left, magenta line); transitions between these two states can be mediated by a Hopf bifurcation.

The simulations performed with our simple model can thus reconcile our experimental observations with the detection of sustained oscillations for different cells (Kellogg and Tay, 2015), but also with the rich variety of damped oscillatory dynamics reported for other cell types in recent works (Cheng et al., 2015; Sung et al., 2014). Other theoretical analyses have identified the amount of NF-κB expressed by a specific cell as the determinant of the Hopf bifurcation and the stability of the oscillator (Mothes et al., 2015). Notably, we used cells with physiological and uniform p65 levels. We also note that the interplay between different feedback loops, including the A20 feedback (Werner et al., 2008) can vary between cell types, leading to different degrees of dampening, as observed experimentally by selective knock-outs (Hoffmann et al., 2002); dissecting the contribution of each of these feedbacks will be the subject of future experimental and theoretical work.

Oscillations in the nuclear concentration of NF-κB lead to different dynamical patterns of gene expression

The oscillatory behaviour of NF-κB consists in periodic relocalization between the nucleus and the cytoplasm, and entails corresponding cycles of NF-κB binding to cognate binding sites in the genome. Thus, NF-κB driven transcription should be pulsatile as well, and in fact we detect cycles in the levels of elongating transcripts in the prototypical early and late NF-κB-controlled genes IkBα and Ccl5. However, whereas mature mRNA oscillates for IκBα, it slowly accumulates for the late gene Ccl5, and this held true for different frequencies of the periodic stimulation. Genome-wide analysis of dynamics, using unsupervised clustering of whole transcriptome profiles, identified six highly homogeneous sets (“clusters”), three containing genes with increasing transcription and three containing genes with decreasing transcription. While the genes with increasing transcription are predicted NF-κB targets, those in the decreasing sets are not.

The train of p65 nuclear translocations is thus decoded by the cell into different transcriptional dynamics: genes in Cluster 1 display a clear and sustained oscillation in transcript level; instead, transcripts for genes in Clusters 2 and 3 accumulate fast or increase slowly and steadily over many cycles of NF-κB oscillations, respectively (see Figure 8C). The result is that each new wave of TNF-α elicits on Cluster 1 genes almost the same responses as the previous waves, and no memory is kept of the past; the expression of genes of Clusters 2 and 3 allows a long-lasting response that integrates over time the responses to previous waves.

All three dynamics of genes with increased expression could be reproduced by our minimal mathematical model of transcription, which incorporates the features of the paradigmatic telegraph model (Suter et al., 2011). Our model also suggests that mRNA degradation is the key parameter that allows a functional and temporal shaping of the NF-κB response, in particular for transcripts in Cluster 1 (Figure 6—figure supplement 6, top-left panel). The computational results fit well with the fast turnover of early transcripts, whose degradation depends on positive regulators like Zfp36 (Rabani et al., 2011). Of note, Zfp36 is in Cluster 1, and oscillates sharply and synchronously with the forcing; this might contribute further to the shaping of peaked responses.

Oscillatory dynamics and functionally related patterns of gene expression: an adaptive viewpoint

Overall, our most remarkable finding is that NF-κB oscillations drive the expression of distinct sets of genes whose expression dynamics and functions correlate: the transcripts that oscillate encode mostly cytokines and cytokine receptors, the ones that rise fast and decrease slowly encode proteins involved in immunity, and the ones that rise steadily encode proteins involved in the rearrangement of the extracellular matrix (Figure 8C).

While steady-rise and rise/decline programs of gene expression can be operated by non-oscillatory transcription factors, oscillating transcription programs can be best achieved with oscillating transcription factors. We note that the transcripts that oscillate are involved in the decision of cells of whether to move or not, which has to be refreshed repeatedly over time. The oscillatory dynamics of NF-κB allows time to be segmented in units of 90 min cycles, with the option to resume chemokine-directed migration after each single cycle. As seen from this perspective, the operation of the NF-κB system as a damped, fast-detuning oscillator implies that in the presence of constant stimulation each motile cell will decide for itself, without synchrony to nearby cells. In contrast, in a situation where stimuli change over time beyond a certain threshold – at least three-fold variation in amplitude – cells would be broadly coordinated. Uncoordinated movement may maximize the volume randomly patrolled by cells, while coordinated movement may be optimal to reach a specific target location.

We thus speculate that NF-κB oscillations may have adaptive value in achieving oscillatory expression of genes involved in movement and direction, and the cells’ decision between uncoordinated or collective action. In conclusion, we suggest that the oscillatory dynamics of NF-κB (and other transcription factors) is a means of segmenting time to provide opportunity windows for decision.

Materials and methods

Please note that this section contains both technical details and supplemental information explaining the rationale behind our experimental approaches. We also refer here to Figure supplements that are also discussed in the main text.

Cell lines and cell culture

GFP-p65 knock-in fibroblasts were provided by M. Pasparakis (details in De Lorenzi et al (2009); Sung et al (2009)) and cultured in phenol-red free DMEM, with 10% FCS, 50 µM β-mercaptoethanol, 1x L-glutamine, 1x sodium pyruvate, 1x non-essential amino acids, and 1x pen/strep in standard tissue culture plastic.

GFP-p65 fibroblast cultures were started from original aliquots frozen upon arrival from Pasparakis' lab and tested every week to exclude mycoplasma contamination using approved kits. Since there are no approved STR profiling protocols for mouse cells ((Yu et al., 2015) and http://www.atcc.org/Global/FAQs/0/A/STR%20Testing%20Service%20-%20STR%20on%20mouse%20cell%20lines.aspx), we carefully tested our cells for the presence of non-green cells (not expressing GFP-p65) that could represent a cross-contamination of the original culture.

For imaging experiments, cells were plated one day before the experiment in CellASIC™ ONIX M04S-03 Microfluidic Plates at low density to avoid confluence on the day of the experiment (e.g. Figure 1—figure supplement 1A). These plates consist of chambers for cell culture connected through microfluidic channels to a series of reservoirs containing media with selected concentrations of stimuli that can be flown through the chambers. Of note, to avoid cell stress or toxicity, the microfluidic plates are primed with 10%FCS in DMEM for 2–4 hr before cell plating.

Before imaging, DMEM-0.1% FCS medium containing 50 ng/ml of the nuclear dye Hoechst33342 was replaced in the microfluidic chambers 3 hr prior to the experiment using the microfluidic platform, see details on the use of the CellASIC™ ONIX below. Mouse recombinant TNF-α (R&D; Systems) was diluted in DMEM-0.1%FCS as specified in the text and added in the plate reservoirs. It is relevant to note that TNF-α activity is maintained after 12 hr incubation in the plate reservoirs at 37°C, indicating that TNF-α degradation is negligible when it is not in contact with the cells (Figure 1—figure supplement 1D).

Microfluidics

The CellASIC® ONIX Microfluidic Platform allows to apply sharp stimuli by quickly replacing the medium in the cell chambers. The flexible proprietary software allows the delivery of medium containing different concentrations of TNF-α for specific time sequences for more than 10 hr while cell are imaged. The protocols are available upon request. The medium flows in the channels around the microfluidic chambers and diffuses through the perfusion barrier, minimizing the undesirable effect of shear stress (Babini et al., 2015) (see Figure 1—figure supplement 1A,B). The small volume of medium in the chambers (less than 1 µl) is replaced fast even at low pressures, as confirmed by the sharp oscillatory profiles obtained even for high frequency stimuli. In spite of the low pressure applied (1 psi), the flow rates of 10 µl/hr (manufacturer's website) are in the same order of magnitude as the 200 nl/min used in Kellogg and Tay (2015). We also generated “sawtooth-like” profiles as in Kellogg and Tay (2015) by periodically replacing the medium with a 15 min pulse and stopping the flow. However, the cell-to-medium ratio in our culture chambers might be different and relevant for TNF clearance due to cell-linked decay needed for the sawtooth profile.

Live imaging and detection of NF-κB dynamics

Live cell imaging of GFP-p65 knock-in MEFs was performed using a Leica TCS SP5 confocal microscope with an incubation system where cells were stably maintained at 37°C in 5% CO2. Time-lapse images were acquired at 6 min intervals for more than 15 hr. We used a low magnification objective (20x, 0.5 NA) and an open pinhole (Airy 3), ensuring that the image depth (10.7 µm) contains the thickness of the whole cell so that images are a record of the total cell fluorescence. GFP-p65 is imaged with the 488 nm Argon laser (GFP channel) while Hoechst 33342 stained nuclei are imaged with the low energy 405 nm UV diode laser at 5% of its maximum intensity (HOE channel). Images were acquired as 16 bit, 1024×1024, TIFF files. Experiment replicates were acquired on different days starting from different batches of frozen cells (samples provided as Videos 15).

Nuclei were stained with Hoechst 33342 for imaging, segmentation and tracking (Zambrano et al., 2014a) without any interference with the natural signalling dynamics (see below for a description of the controls performed).

Safety assessment for Hoechst staining

Stress from photo-damage during live-cell imaging can hamper cell biology studies (Cole, 2014), while nuclear dyes can interfere with the natural signaling dynamics (Ge et al., 2013; Martin et al., 2005). Hence different controls were performed to exclude distortions of the natural signalling in our imaging conditions.

Apoptosis and cell death

We noticed that cells under a constant flow of fresh medium do not display apoptotic events until very late in the experiments (10–12 hr; Video 2). Importantly, unstimulated cells are viable and undergo mitosis despite the presence of low serum concentration (0.1% FCS). In contrast, cells receiving a constant flow of TNF in the absence of Hoechst and UV imaging undergo apoptosis (Video 1), with a frequency that is proportional to the TNF concentration (Figure 1—figure supplement 7). The imaging fields contain detectable apoptotic cells after 3, between 6 and 9 or after 12 hr of stimulation with 10, 1 or 0.1 ng/ml of TNF, respectively. This suggests that cell death observed in Video 1 was due to the continuous flow of 10 ng/ml TNF-α and not due to the imaging conditions.

Photo-damage

The cells were either imaged for 15 hr in the presence of Hoechst or left without Hoechst and not imaged; we then immunostained cells to detect thymine dimers formation. Immunostaining intensities (Figure 1—figure supplement 6) are extremely low in cells exposed to Hoechst and imaging, especially when compared to cells exposed to UVC radiation. This suggests that imaging does not produce thymine dimers or that thymine dimers are repaired fast and do not accumulate in the cells.

To further exclude the possibility of DNA damage, we checked for the presence of gammaH2AX, which is a marker for early DNA repair after formation of DNA double strand breaks (DSB) (Turinetto and Giachino, 2015; Yang et al., 2015) and late DNA repair after the formation of UV-induced thymine dimers (Oh et al., 2011; Staszewski et al., 2008).

Again, the signal for gammaH2AX (Figure 1—figure supplement 7) is extremely low in imaged cells, when compared to the signal from cells treated with doxorubicin, which induces DSBs. Of note, the phosphorylated form of H2AX is found at DNA polymerase stalled forks (Turinetto and Giachino, 2015) as marker of replicative stress. This observation may explain the minimal signal present in some of the imaged cells. All together these results indicate negligible DNA damage in the conditions considered.

Interference with NF-κB dynamics

We performed a manual tracking of cells that received a constant flow of TNF-α but were not exposed to Hoechst nor to UV light (“unstained/non-imaged cells”, Figure 1—figure supplement 9). Manual segmentation of nuclei is simple when nuclei are either empty or full of p65, but becomes unreliable when concentrations in nucleus and cytoplasm are similar. Despite this caveat, we found a good qualitative agreement of dynamics in unstained/non-imaged cells with the heterogeneous but damped oscillatory dynamics observed under Hoechst/UV imaging conditions, as shown in Figure 1—figure supplement 9. Therefore, the p65 dynamics described in our manuscript represents the biological effect of stimulation with TNF-α, and is essentially undistorted by imaging.

Immunobloting and immunostaining

Staining was performed according to manufacturers’ instructions. For WB, we used an anti p65 Rab (dil 1:1000, #sc-372 C20, Santa Cruz; Figure 1—figure supplement 10). For immunostaining we used anti gammaH2AX Mab (dil. 1:500; #05-636, Millipore; Figure 1—figure supplement 7) and anti thymine dimers TDM-2 (dil. 1:2000, Cosmo Bio, Japan (Komatsu et al., 1997); Figure 1—figure supplement 6). UV-photodamage was induced with 254 nm UV irradiation using an UVC 500 Crosslinker, Amersham. Panels in Figure 1—figure supplements 6 and 7 are representative of 3 independent experiments.

Quantification of NF-κB dynamics

The following short description summarises the whole process in few sentences. Each passage will be then thoroughly described in the next paragraphs.

The dynamics were quantified by computing the nuclear to cytoplasmic ratio of the intensities (NCI) of single cells, which is a measure robust against distortions and reflects faithfully the oscillations in the nuclear concentration of NF-κB oscillatory dynamics. This holds true provided that the total amount of p65 is constant, as we show for our cells under stimulation (see Figure 1—figure supplement 9), and is also assumed by mathematical models (see Hoffmann et al (2002) and the more recent in Kellogg and Tay (2015)). Interestingly, p65 is not constant in macrophages stimulated with LPS (Sung et al., 2014). The analysis of time series relies on the detection of significant peaks, performed following our previously discussed procedure (Zambrano et al., 2014a) that allows to distinguish significant peaks from noisy peaks. Once peaks are detected we assign them a phase: 2π for the “maxima” and π for the “minima” between peaks.

Selection of the quantifier of NF-κB dynamics

In order to assess the oscillations for a larger number of cells, we optimized our software (Zambrano et al., 2014a) to calculate the nuclear to cytoplasmic ratio of the intensity (NCI) of NF-κB signal for hundreds of cells; NCI is a quantifier that has been used by other groups (Ashall et al., 2009; Nelson et al., 2004). As we argue below, thanks to the fact that the total amount of NF-κB is constant (Figure 1—figure supplement 10), NCI depends univocally and monotonically on the nuclear amount of NF-κB. More importantly, due to the fact that it is a ratio of intensities, it is robust and independent of slight changes in the focus and in the laser intensity, among other possible experimental distortions, which are observed in our setup. The same rationale brought us to use ratios of intensities in our previous works (Sung et al., 2009; Zambrano et al., 2014b). Importantly, these distortions imply that the background-adjusted mean nuclear intensity used in other studies (Kellogg and Tay, 2015; Lee et al., 2014; Sung et al., 2014), although advantageous for other reasons (it only requires the segmentation of the nuclei and estimation of the background) would not be appropriate in our imaging experimental setup.

We provide below a more detailed argumentation of these ideas.

Following the notation of Zambrano et al (2014b) we have that the intensity measured in pixel p of the image at time t in a time lapse experiment can be described as:

(Q1) I(p,t)=A(p,t)P(p,t)+Bp,t

where P(p,t) is the amount of NF-κB in pixel p, A(p,t) is the amplification coefficient between the protein brightness and the amount of protein and B(p,t) corresponds to the background. In our experiments with time-varying intensities, it is clear that both A(p,t) and B(p,t) vary in time, and presumably also in space, due to laser variations and/or slight variations in the illumination uniformity.

Our quantifier NCI, that we estimate using our software, is the ratio of the background corrected average intensities in the nucleus and the cytoplasm I(t)nuc and I(t)cyto respectively, and can be defined using this notation as:

(Q2) NCI(t)I(t)nucI(t)cyto=ScytoSnucleuspnucleusI(p,t)-B(p,t)qcytoplasmI(q,t)-B(q,t)=ScytoSnucleuspnucleusA(p,t)P(p,t)qcytoplasmA(q,t)P(q,t)

so we have that, if we consider that A(p,t) is approximately constant and equal to certain A(t) in the area occupied by the cell (as it is in our images):

(Q3) NCI(t)ScytoSnucleusA(t)pnucleusP(p,t)A(t)pcytoP(q,t)=ScytoSnucleus(NFκBnuc(t))(NFκBcyto(t))

where Scyto and Snuc are the areas occupied by cytoplasm and the nucleus, respectively. Our quantifier is hence a good indicator of the ratio of the amount of NF-κB in the nucleus and in the cytoplasm. However the numerator and denominator of this ratio can fluctuate due to biological reasons and this might blur the existence of oscillations in the nuclear amount of NF-κB. But this is not the case, due to the fact that the total amount of NF-κB is constant for our stimulated fibroblasts. This has been assumed in different mathematical models present in the literature, from the seminal paper (Hoffmann et al., 2002) to more recent papers (Kellogg and Tay, 2015). Interestingly, though, this view has been challenged by the recent work of Sung et al. (2014), which shows that in macrophages under LPS stimulation p65 is regulated by a positive feedback. Thus, to confirm that our assumption is reasonable, we quantified p65 by western blotting for our cells under 10 ng/ml TNF-α at several timepoints from 1 to 8 hr. The results are shown in (Figure 1—figure supplement 10), showing no significant change in the p65 levels upon stimulation. Hence, we consider our assumption valid. This implies that:

(Q4) NCI(t)=ScytoSnucleus(NFκBnuc(t))(NFκBTOT)(NFκBnuc(t))

hence we can see that NCI(t) depends monotonously on (NFκBnuc(t)), see Figure 1—figure supplement 9B. For this reason, it is intuitively clear that each local maximum or minimum of (NFκBnuc(t)) leads to a local maximum or minimum of NCI(t). Hence, oscillations in the nuclear amount of NF-κB will be observed also using NCI(t). We can put this mathematically by saying that oscillations in the nuclear concentration of NF-κB occur at times t for which the condition (NFκBnuc(t))'=0. Similarly, oscillations in NCI will depend on the value of the derivative of NCI, that is:

(Q5) NCI'(t)((NFκBTOT)(NFκBnuc(t))+1)((NFκBTOT)(NFκBnuc(t)))2(NFκBnuc(t))'

so from the above formula it is easy to see that:

(Q6) NCI'(t)=0(NFκBnuc(t))'=0

Overall, then, imaging, mathematical and biochemical arguments confirm that computing NCI(t) is an adequate way to quantify oscillatory dynamics.

Concerning the background-corrected average intensity, using the above notation it would be calculated as:

(Q7) I(t)nucleus=1SnucleuspnucleusI(p,t)B(p,t)=1SnucleuspnucleusA(p,t)P(p,t)

thus I(t)nucleuspnucleusP(p,t)=NFκBnuc(t) for all t only if A(p,t) were constant in all time frames. As argued previously, this is not the case in our setting and that’s the reason why we preferred to use NCI.

Finally, variations in the areas of the nucleus and the cytoplasm might introduce small distortions, although they typically vary seldom and slowly. Notice though that our software discards cells for which the areas change abruptly, which typically is an indicator of imminent mitosis or cell death. We cannot totally exclude variations of p65 at single cell level, provided that we measured it using a population assay; however, p65 turnover is known to be long, so its slow variation would contribute with a small term in the derivative and thus would just slightly distort the times for which peaks in NCI appear compared to those of the nuclear concentration. To conclude, we think that NCI would still capture the peaks and hence would be able to assess the cell's oscillatory state, which is our aim here.

The data provided confirm the validity of NCI as a quantifier. First, we have that NCI(t) remains reasonably constant – except for some spontaneous oscillations, as previously reported (Zambrano et al., 2014b) for cells that are not stimulated, as shown in Figure 1—figure supplement 9A. This is remarkable because the movie from which these series were obtained present considerable variation in the image intensity (see the Video 2). Along the same lines we do not find upwards or downwards average trends in our data of NCI time series for different conditions discussed in the main figures and in the figure supplements, which indicates further that the time series are properly normalized and reflect well the dynamics. As an additional confirmation, we have plotted together the average nuclear intensity and the NCI values obtained from manual segmentation in Figure 1—figure supplement 9C, that show the same kind of qualitative behavior as NCI but with different trends, as expected. A last indicator of the sensitivity of the quantifier comes from its ability to detect even small oscillatory peaks, see Figure 3 and Figure 3—figure supplement 1 as an example. Overall, then, NCI is a faithful measure of the oscillatory state of our cells.

Automated analysis of time-lapse experiments data

The major advantage of considering NCI with respect to the nuclear to total ratio used in Zambrano et al (2014a) is that it does not require a perfect segmentation of the cytoplasm, which might be complicated when cells touch each other. We developed a software that uses this quantifier and is thus able to multiply by a factor of 2 the number of cells tracked, and by a factor of 1.5 the average tracking time with respect to the one described in Zambrano et al (2014a).

The software used to calculate NCI is provided as source code and works as follows: for each time-lapse experiment, we have N frames images in the HOE channel and in the GFP channel, respectively. Nuclei were segmented and used for cell tracking following the procedure described in Zambrano et al (2014a). In order to estimate the average cytoplasmic intensity, the background was computed by taking a square area centred on the cell nucleus, dividing it in tiles and using the one with the smallest average intensity in the GFP channel to estimate the background intensity (this procedure gives values compatible with the values obtained using a clustering-based algorithm). Points belonging to the cytoplasm are those around the nucleus in a size window equal to 1.5 times the size of the nucleus. The average cytoplasmic intensity is computed from the intensity in the GFP channel of the pixels from this “ring”. An example is shown in Figure 1—figure supplement 1D.

Dividing or apoptotic cells were identified assessing their geometrical features (abrupt changes in size of the nucleus and of the “cytoplasmic ring”) and discarded automatically.

Time series analysis

To analyse the resulting NCI time series, we adapted our peak detection algorithm (Zambrano et al., 2014a). We consider peaks as a sequence of a local minimum, local maximum and local minimum. The value of the peak θ is defined by the height of the peak (from the highest minimum to the maximum). A peak defined by only three consecutive time points is considered a noise peak. We plot the distribution of the noise peak values obtained from our unperturbed and chronic stimulation time series (Figure 1—figure supplement 2A), for which we observe our previously reported spontaneous activations (Zambrano et al., 2014a). We find that these noisy peaks have a low value, and hence by considering as significant those with a value over θ>0.15, we find a reasonably good compromise between the need to ignore noise peaks and the need to detect small peaks of valuable dynamical information. This was tested in a number of time series, and an example of time series with the significant peaks detected using this threshold is shown in Figure 1—figure supplement 2B. Calculations of magnitudes inferred from peaks take advantage of time series from cells that were tracked for at least 7 hr.

Following (Mondragon-Palomino et al., 2011) we used the peaks to assign a phase value: peak maxima are assigned a phase 0 (2π) and a phase π is assigned to peak minima. This was of particular interest to obtain an assessment of the oscillatory modes present in our system, provided that peaks can be very heterogeneous and plotting the peak height in colour-plots might make it difficult to appreciate the smaller ones and hence the possible resonant oscillatory patterns. An example of this transformation is shown in Figure 1—figure supplement 2C, which is the phase derived from the peaks obtained of the time series displayed in Figure 1—figure supplement 2B. As in Mondragon-Palomino et al (2011) we use the phase difference between the time series and the forcing to quantify the degree of synchrony. The phase difference Δφ was calculated from the timing ΔT between the beginning of each forcing cycle and the closest significant peak, as Δφ=2πΔT/Tf, where Tf is the forcing period. The entropy of the distribution of the phase is calculated as:

S=kpklogpk

where pk is the probability of the kth bin (we use eight bins for all the conditions considered). The synchrony intensity is then computed as:

η=1S/Smax

where Smax is the value obtained for equal values of pk.

Quantitative real-time PCR

Total RNA was isolated using the IllustraRNAspin Mini kit (GE Healthcare), and complementary DNA (cDNA) was obtained by retro-transcription with Random Hexamers and SuperScript II Reverse transcriptase (Invitrogen) following the manufacturers’ instructions. Primers for the detection of mature transcripts were designed in adjacent exons and spanned the intervening intron; primers for nascent transcripts were located across an exon-intron boundary. Primer sequences are listed in Supplementary file 1.

Quantitative real-time PCR was performed with SYBR Green I protocol using the LightCycler480 (Roche). The results are shown as averages of three technical replicates. Analysis was performed with the –ΔCt method corrected for primer efficiencies (Vandesompele et al., 2002) and normalised with two reference genes (Actb and Rplp1) (Nordgard et al., 2006). Experiments were repeated twice.

It is important to emphasize that we did not stain our cells nor did we illuminate them with our UV laser in our transcription experiments, neither in the RT-PCR assays nor in the Microarray experiments described below. The controls described previously show that the nuclear labelling and the imaging do not interfere with NF-κB signalling dynamics, so we expected the same for NF-κB-driven transcription. To further confirm this, we compared RT-PCR quantifications of transcription in cells stimulated with TNF and cultured in imaging medium (0.1% FCS+Hoechst) or standard MEF medium (10% FCS). The results for nascent and mature IκBα and Ccl5 transcripts reported in Figure 5—figure supplement 2 show no difference in the transcriptional response in the two conditions.

Microarray experiments

RNA samples were extracted using the IllustraRNAspin Mini kit (GE Healthcare). Following extraction, RIN (RNA Integrity Number) was >9 (BioAnalyser, Agilent RNA Nano Kit). RNA samples (500 ng) were reverse transcribed with the IlluminaTotalPrep RNA Amplification Kit (Ambion) and copy RNA (cRNA) was generated with 14 hr in vitro transcription reactions and checked at the BioAnalyser. Washing, staining, and hybridization were performed according to standard Illumina protocols. cRNA samples were then hybridized to IlluminaBeadChip Array MouseRefSeq-8 v2. BeadChips were scanned with BeadArray™ Reader in channel 2. The data have been uploaded on the Dryad Digital Repository (Zambrano et al., 2016). Experiments were repeated twice.

Bioinformatic analysis of microarray experiments

Genome Studio’s bead summary probe level data – not normalized and not background corrected – were analyzed using Bioconductor. We performed quality assessment by plotting the intensities of regular and control probes. Sample intensities were quantile normalized and filtered for expression and probe quality using beadarray R package: only probes with detection P-value <0.05 in at least one condition and whose categories is defined “Perfect” or “Good” were kept. Probes were labelled as deregulated if their absolute log2 fold change relative to the 0 time point were greater than 1.

Clustering was performed applying a soft clustering of deregulated probes with the Mfuzz R package, which is suggested for microarray time course-data (Kumar and Futschik, 2007). The optimal number of clusters was assessed with the d.min function. To focus on the shape of the gene expression profiles, we standardised the gene expression profiles by taking the usual base-two logarithm and normalizing to obtain mean 0 and standard deviation 1. The algorithm then groups genes based on the Euclidean distance between profiles and the c-means objective function, which is a weighted square error function. Each gene is assigned a membership value between 0 and 1 for each cluster. Hence, genes can be assigned to different clusters in a gradual manner. The parameters m defines the degree of “fuzzification”. It is defined for real values greater than 1 and the bigger it is the more fuzzy the membership values of the cluster. We used an m value estimated by the “mestimate” function of 1.5. Membership values indicate the similarity of vectors to each other defining a cluster cores. To extract list of genes belonging to the cluster cores, we used the “acore” function taking from each clusters all genes with a membership value of at least 0.5.

Pathway analysis of genes contained in each cluster was performed with the ClusterProfiling R packages using an adjusted p-value cut-off for enrichment <0.2 based on the hypergeometric distribution (Yu, 2015). We used for our analyses the Reactome database. In order to simplify the resulting output we plotted only the top five categories of the second and third level (Yu, 2015). Conversely, overlaps were calculated based on all categories of the second and third levels. The overlap coefficient (or Szymkiewicz-Simpson coefficient) between gene sets X and Y is given by:

overlap (X,Y)=|X  Y|min(|X|, |Y|)

For statistical significance we performed a Fisher's Exact Test for evaluating if the resulting clusters where enriched for NF-κB's targets. The resulting p-value has been converted in a significance measure (-Log10(p value)). We used both a list taken from Li et al (2014) and from Thomas Gilmore’s website, Boston University (http://www.bu.edu/nf-kb/gene-resources/target-genes/).

Mathematical modelling

We propose here a mathematical model based on the one discussed in Zambrano et al (2014b) that adds the layer of regulation by considering the negative feedback provided by the protein A20 that blocks IKK activation upon stimulation (Figure 1—figure supplement 4). For the sake of completeness, we describe briefly below the basic process that we considered, together with the biochemical rates involved (values are given in Supplementary file 2) and a summary of our normalization procedure. Additional details on the motivations underlying certain selection of biochemical reactions and variables of interest can be found in that paper. Being a simple model, it provides important qualitative insights on the possible variety of single-cell dynamics, but we also use it to provide quantitative fittings of the average population dynamics and transcription.

NF-κB activation and IκBα feedback

The biochemical reactions described here are essentially the basic ones given in the simple model given in Zambrano et al (2014b). The basic simplification of our model is to consider that free NF-κB is nuclear, while the complex with the inhibitor IκBα is cytoplasmic. The formation of the complex is given by (NF-κB:IκBα)

NFκB+IκBαA(NFκB:IκBα)

that can also spontaneously dissociate

(NFκB:IκBα)dNFκB+IκBα

An external signal can lead to the appeareance of IKKa that can free NF-κB by degrading of IκBα in the complex

(NFκB:IκBα)+IKKaPNFκB+IKKa

and in its free form

IκBα+IKKaκPIKKa.

The negative feedback loop is enabled by the fact that the gene encoding IκBα, Gα, can be activated by NF-κB

Gα,OFF+NFκBKon,IGα,ON+NFκB

while the inactivation is modulated by

Gα,ON+IκBαKoff,IGα,OFF+IκBα.

Notice that we are assuming here that IκBα plays a role as transcriptional repressor, as suggested by experimental results (Arenzana-Seisdedos et al., 1995) and used in some mathematical models (Kellogg and Tay, 2015; Lipniacki et al., 2007; Tay et al., 2010), but not others (Ashall et al., 2009; Sung et al., 2014). We show later that this does not have a strong impact on the fittings, but to facilitate the correspondence with our previous work and in line with the mentioned modelling approach, we opt to keep this process for all the gene inactivations considered in most of our explorations. We only briefly explore the effect of setting the koff=0 in Figure 6—figure supplement 7, results are detailed in the text.

We also consider that the gene can be basally activated

Gα,OFFkon0,IGα,ON

and inactivated

Gα,ONkoff0,IGα,OFF

Transcription (mRNA production) is given by

Gα,ONKRIGα,ON+IκBαRNA

while the RNA degradation is given by

IκBαRNAdRI.

IκBα translation is given by

IκBαRNAKIIκBα

while its spontaneous degradation is given by

IκBαdI,

a degradation that is also possible while it is forming the complex

(NFκB:IκBα)γdINFκB

IKK activation and A20 feedback

The protein A20 is known to provide a negative feedback by modulating IKK activation process (Ashall et al., 2009). We summarize this in our new mathematical model using the following processes in which we just model the evolution of the active form, IKKa.

Given a external time dependent TNF-α variation, TNF(t), characterized by certain values of the alternating doses, Dand D2 in times Tand T2, we will assume that the appearance of the active IKK, IKKa, can be summarized as:

TNF(t)K(A20)IKKa.

In a situation of constant flow of TNF-α the value of TNF(t) would be constant and equal to a dose D.

The biochemical rate of this equation, K(A20), is the only rate that we consider as variable. In doing so, we aim to summarize in just one biochemical reaction the fact that A20 contributes to block IKK activation (Ashall et al., 2009) instead of modelling the whole IKK activation module, so

K(A20)=KS1+A20A200n

IKK gets spontaneously inactivated as

IKKadK.

The gene encoding A20 is regulated by the same activation and inactivation processes as IκBα so

GA20OFF+NFκBKon,A20GA20,ON+NFκB
GA20,ON+IκBαKoff,A20GA20,OFF+IκBα
GA20,OFFKon0,A20GA20,ON
GA20,ONKoff0,A20GA20,OFF

Transcription is given by

GA20,ONKR,A20GA2,ON+A20RNA

while the RNA degradation is given by

A20RNAdR,A

finally A20 is translated as

A20RNAKAA20

and it is spontaneously degraded as

A20dA

Transcription of an NF-κB controlled gene

As for the genes encoding for A20 and IκBα, we consider that for any NF-κB controlled gene transcription is regulated by the processes

GOFF+NFκBKon,GGON+NFκB
GON+IκBαKoff,GGOFF+IκBα
GOFFKon0,GGON
GONKoff0,GGOFF

while the RNA produced from the gene is provided by

GONKRGON+A20RNA

and its degradation is given by

GRNAdR,G

Model equations, normalization and parameter selection and uncertainty

The equations of the model are derived using mass-action kinetics and performing a normalization in much the same way as in Zambrano et al (2014b). We describe below the process, using the same symbols to represent biochemical species and their copy number.

First, as in previous existing models and as suggested by our experiments, the total amount of NF-κB, free plus bound to the inhibitor, remains constant and equal to NF-κB0 so we define N=NF-κB/ NF-κB0. We also normalize the amount of inhibitor as I= NF-κB/ NF-κBand the amount of A20 as A= A20/ A200 and K=IKK/IKK0, where A20and IKK0 are amounts of reference in such a way that N, I and A are adimensional variables of the order of 1. On the other hand, for any gene we have that the maximum number of on states is G0=2 (the two alleles), so GI, and GA are a number between 0 and 1 representing the fraction of active genes encoding for IκBα and for A20, respectively – it is a continuous approximation to a discrete variable. Finally, notice that using our notation the maximum number of active genes is G0 so it can be proved that the maximum asymptotic value of the amount of RNA of the biochemical species X is of the form KR,X·G0/dR,X. Hence, we define the variables RI, and Ras proportional to the amount of mature RNA of IκBα and A20, respectively, via the relations

RI=IκBα RNA ·dR,I/(KRI·G0) and RA=A20RNA ·dR,A20/(KR,A20·G0).

Overall, these normalizations allow one to pass from biochemical species with copy numbers of different orders of magnitudes to adimensionalized variables that are of the order of 1.

By using them, and renaming combinations of constants (in such a way that lowercase parameters are obtained from uppercase biochemical reaction rates), we have the following model of ordinary differential equations describing the dynamics of our system:

(1) dKdt=dK·K+11+(A)nS(t)
(2) dNdt=d·(1N)+γ·dI(1N)+p·K·(1N)a·I·N
(3) dGIdt=(konIN+kon,0,I)(1GI)(koffII+koff,0I)GI
(4) dRIdt=dRI(GIRI)
(5) dIdt=d·(1N)κ·p·K·Ia·I·N+kI·kRIdI·I
(6) dGAdt=(kon,AN+kon,0,A)(1GA)(koff,AI+koff,0,A)GA
(7) dRAdt=dRA(GARA)
(8) dAdt=kA·RAdA·A

The parameters used as starting point in our explorations and the fittings can be found in the source code folder for mathematical modelling, while the original biochemical rates from which they were derived can be found in Supplementary file 2. Normalised parameter values are provided in the source codes. Many of them have the same values as the ones provided in Zambrano et al (2014b). Additional parameters were manually fitted or extracted from (Tay et al., 2010), as we did previously. Notice that the degradation rates of the kinase, the RNAs and the proteins appear as parameters in the equations above, unchanged. The remaining parameters are related to the original biochemical rates as p=P·IKK0,a=A·NF-κB0, kI=KI·KR,I·G0/(dR,I·NF-κB0), kA=KA·KR,A·G0/(dR,A·A200), kon,I=Kon,I· NF-κB0, koff,I=Koff,I· NF-κB0, kon,A=Kon,A· NF-κB0, koff,A=Koff,A· NF-κB0. S(t) is the normalized signal, in such a way that a constant TNF = 10 ng/ml is equivalent to S(t)=2 h-1 in the system of equations. For lower doses we use lower values of S.

Notice that the adimensionalization leads to a number of parameters smaller than the total number of biochemical rates and constants of the original system. Equations 2–5 are identical – except for the inclusion of spontaneous gene activation-inactivation – to the ones proposed in our previous minimal model of the regulatory network (Zambrano et al., 2014b). Equation 1 models the effect of A20 as inhibitor of the kinase activation. This second negative feedback is completed with Equations 6–8, that explicit NF-κB control of A20 expression.

As shown below, we use this model both to fit the observed dynamics and for numerical exploration of the dynamics observed under different stimuli. For simplicity, we call PNF-κB to the set of parameters of Equations 1–8 and PS to the parameters describing the parameters of the external signal, that can either be constant or a time varying rectangular forcing signal as shown in Figure 1A. Following what we did in previous numerical explorations (Zambrano et al., 2014b) we associate to each parameter an uncertainty degree D, so each parameter can be randomly varied by multiplying it by a factor in the interval [10-D, 10D] with D smaller or equal than 1. In other words, parameters with higher uncertainty degree can be varied up to one order of magnitude above or below their selected initial value, which is itself a moderate variation. Those degrees are specified in Supplementary file 2: note that we choose higher values of D for manually fitted parameters and lower for those from the literature or from our own measurements. We also consider a high value of D for parameters that summarize many different processes or for those in which our model is less detailed, as in the IKK activation – A20 regulatory module.

Finally, the equations for a gene under the control of NF-κB are:

(9) dGdt=(kon,GN+kon,0,G)(1G)(koff,GI+koff,0,G)G
(10) dRdt=dR,G(GR)

where in this case R=G RNA ·dR,G/(KR·G0), kon,G=Kon,G· NF-κB0, koff,G=Koff,G· NF-κB0. The parameters in (9) and (10) are used to model and fit different expression profiles are denoted as PG. Notice that to model the process of gene expression one needs to set PS, PNF-κB and PG. In this context, two different genes expressed under the same external stimuli conditions would share the parameters PS and PNF-κB but different values of PG. This idea is used in our fittings of dynamics and transcription in this work.

Classification of the different dynamical responses for constant stimulus

We have observed experimentally that even for constant stimulation different dynamics are possible, with trajectories showing different degrees of dampening and sometimes even irregular and non-oscillating profiles (Figure 1—figure supplement 2D) as pointed out in previous works (Sung et al., 2009; Zambrano et al., 2014b). Our relatively simple model illustrates why this situation might arise and how a wider variety of dynamics might arise compared to the regular oscillations documented in the literature for more complicated models (e.g. (Ashall et al., 2009; Tay et al., 2010)).

To illustrate this, we analyzed the dynamics for a constant external stimulus (S(t)=2 h-1) and varying simultaneously each parameter of the parameter set PNF-κB according to their uncertainty degree D. We used them to generate trajectories and selected those giving rise to a response, i.e. those for which N(t) is bigger than 0.4 in the first 3 hr and whose average values in the last hours was smaller or equal than N=0.4. To characterized in a systematic way their dynamics, we located the fixed point of Equations 1–8 and calculated the eigenvalues {λi} (i=1,...8) of the Jacobian. In Figure 1—figure supplement 5A we plot the real and imaginary parts of each set of eigenvalues. We represent with a red dot the eigenvalues belonging to a set in which the real part of at least one of them is bigger than zero. This means that the fixed point for those parameters is unstable and hence trajectories converge to a stable limit cycle around it. In other words, for those parameter combinations sustained oscillations arise. The percentage of parameter sets giving sustained oscillations is below 10% (Figure 1—figure supplement 5B, right panel), which make us conjecture that parameters giving rise to sustained oscillations are not the norm.

We also plot in Figure 1—figure supplement 5C the parameter values giving rise to oscillating and non oscillating (but responding) trajectories. We note that the intervals are very similar except for parameters n and dA involved in the A20 negative feedback, which indicates that this feedback is critical in order to obtain a sustained or a damped oscillatory response. In previous works it was shown experimentally and through mathematical models that the interplay between the IκB and A20 regulatory modules regulate different phases of NF-κB response (Werner et al., 2008). Our numerical results suggest that it also plays a key role in the type of dynamics observed (oscillatory versus damped). Overall, this analysis further hints to the fact that it is the precise combination of parameters, rather than the precise single values, what determines the type of dynamics observed and might account for the different dynamics observed with respect to other groups (Ashall et al., 2009; Cheng et al., 2015; Kellogg and Tay, 2015; Lee et al., 2009; Sung et al., 2009; Tay et al., 2010).

Finally, we have to notice that being the system given by Equations 1–8 a dynamical system with real variables, the imaginary eigenvalues come in complex conjugate pairs. The distribution of parameter combinations with 2, 4 and 6 imaginary eigenvalues are shown in Figure 1—figure supplement 5B (left panel). Notice that an oscillatory frequency can be associated to each couple of complex conjugate pairs (Goldstein et al., 2001), so having more than one of such couples means that the dynamics will combine different frequencies.

To further illustrate all this, in Figure 1—figure supplement 5D we present examples of the variety of trajectories found. Some of them are oscillating (red), showing both smooth and spiky oscillatory peaks. Many others are damped and in some cases the concurrence of two oscillating timescales, fast and slow, is evident. Some others present clearly non-oscillating profiles. Overall, our numerical exploration of our relatively simple models shows that different dynamics are possible in presence of constant stimulation, as observed in the experiments of this and previous works (Sung et al., 2009; Zambrano et al., 2014a).

The model simulation routine was implemented in C code language and compiled using gcc. Routines for parameter variation and time series analysis and stability analysis (using the function fsolve of the Nonlinear Equations package) where written using GNU Octave.

Parameter fitting

Routines were written in GNU octave to fit the model to the data by combining Markov Chain Monte Carlo for initial exploration of the parameter space and a Levenberg–Marquardt algorithm. Parameters of Equations 1–8 were varied within the specified uncertainty degree. The goal is to minimize the distance between a given goal signal X={X(tk)} for k=1,...,N and the theoretical values x={x(tk)} obtained from the model. We define the distance of a given fitting to the data as:

dX,x=1NkX(tk)xtkmax(X(tk),x(tk))

When several data are fitted simultaneously, the total distance between model predictions and the data were obtained by summing each distance. Notice that the 1/N factor weighs for difference in the length of the data. For a given computations, this distance gives the average relative error per timepoint of the fitting.

- To fit NCI, theoretical NCI profiles from the model were obtained as:

NCI(t)=ScytoSnucleusN(t)1N(t).

where Scyto/Snuc is the ratio of the nuclear to total area of the cells in our image, which is estimated to be of around 1/3.

For parameter fitting of the gene expression obtained from Quantitative Real-Time and microarray experiments (see below), the fold change levels predicted by the models where computed as R(t)/R(0). For the simultaneous fitting of NCI, IκBα and Cccl5 mRNA shown in Figure 5, common PNF-κB parameters where found while for IκBα and CCL5 mRNA we independently fitted the respective set of parameters PG. For the fitting of PNF-κB the parameters were varied within the limits given by their uncertainty degree.

Fitting of the microarray data was performed by keeping constant the parameters PNF-κB obtained fitting the data of Figure 5 and varying PG (those of (9–10)) for each gene, as we did for IκBα and Cccl5. Examples are shown in Figure 6—figure supplement 6. Figure 6—figure supplement 5 shows that the distance between fitting and real data is better for genes with increased transcription with an average fitting error of 9%.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
    Ultraviolet radiation-mediated damage to cellular DNA
    1. J Cadet
    2. E Sage
    3. T Douki
    (2005)
    Mutation Research/Fundamental and Molecular Mechanisms of Mutagenesis 571:3–17.
    https://doi.org/10.1016/j.mrfmmm.2004.09.012
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
    Live-cell imaging
    1. R Cole
    (2014)
    Cell Adhesion & Migration 8:452–459.
    https://doi.org/10.4161/cam.28348
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
    Classical Mechanics (3rd Edn)
    1. H Goldstein
    2. CP Poole
    3. JL Safko
    (2001)
    Essex: Pearson Education Limited.
  20. 20
  21. 21
    RNA splicing regulates the temporal order of TNF-induced gene expression
    1. S Hao
    2. D Baltimore
    (2013)
    Proceedings of the National Academy of Sciences of the United States of America 110:11934–11939.
    https://doi.org/10.1073/pnas.1309990110
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
    Synchronization
    1. A Pikovsky
    2. M Rosenblum
    3. J Kurths
    (2003)
    A Universal Concept in Nonlinear Sciences, Cambridge: Cambridge University Press.
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72

Decision letter

  1. Suzanne Gaudet
    Reviewing Editor; Dana-Farber Cancer Institute, Harvard University, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your work entitled "NF-κB oscillations translate into functionally related patterns of gene expression" for peer review at eLife. Your submission has been favorably evaluated by Aviv Regev (Senior editor), a Reviewing editor, and two reviewers.

The reviewers have discussed the reviews with one another and the Reviewing editor has drafted this decision to help you prepare a revised submission

Summary:

How transcription factors affect responses to extracellular signals remains a key question in systems biology for which the NF-κB signaling pathway in mammalian cells has been a canonical study system. The manuscript "NF-κB oscillations translate into functionally related patterns…" by Zambrano et al. presents a systematic investigation of single-cell NF-κB signaling dynamics and of population-level transcriptional responses in GFP-RelA knock-in mouse fibroblasts stimulated with controlled TNF pulses. The authors use a microfluidic device to deliver TNF to cells while avoiding the accumulation of secreted proteins. Although, this does not exclude cell-centric feedback regulation of the signaling response to TNF, the minimal ODE model proposed by Zambrano et al. nonetheless accounts well for their experimental observations and the authors conclude that NF-κB behaves as a damped oscillator in these fibroblasts. That said, the reviewers have raised serious concerns about the single-cell data (detailed below) that need to be addressed.

The second half of the study, focused on the gene expression kinetics and functional correlates in cells treated with TNF pulses, leverages the use of oscillatory stimulations to reveal different clusters of genes that share dynamic patterns of expression. This is the most convincing and interesting finding of the manuscript. An important biological insight is that different functional pathways exhibit different expression dynamics; chemokine/chemokine receptors oscillate with TNF inputs, while immune responses integrate over longer timescales. The authors convincingly show that the decay rate of the mRNA of a gene of interest is sufficient to account for the different dynamic regimes of gene regulation.

Essential revisions:

1) The live cell imaging conditions seem to have been very damaging to cell physiology. Stress from photo-damage during live-cell imaging is a very serious issue which has hampered many cell biological studies (Cole R, 2014 Cell Adhesion & Migration 8 (5) 452-259), and it is essential to address this problem. In the data presented by Zambrano et al., the supplementary movies show that the knock-in fibroblasts express very low levels of GFP-RelA making imaging challenging and it is quite evident from these movies that there is significant cell death and stress, and few cell divisions. Photo-stress was likely increased in this study due to the repeated exposure to the UV laser used for Hoechst imaging. A key concern is that with such a high level of stress, the observed NF-κB dynamics may not represent the true signaling dynamics induced solely by TNF. This could explain the authors observe highly damped oscillations, which are different from the sustained oscillations observed and reported in their first paper on the same GFP-RelA knock-in MEFs (Sung MH et al. 2009 PLoS One). In contrast to the manuscript under consideration, no UV or Hoechst were used in the 2009 study and healthy cell divisions were frequently observed throughout the long-term imaging supplementary movies. In light of these important concerns, the inclusion of the single-cell imaging data in the paper requires additional control data:

A) Because nuclear dyes can alter cell cycle and/or cause oxidative stress and mutations (Martin RM, Leonhardt H, and Cardoso MC, 2005 Cytometry; Ge J, Wood DK, et al. 2013 Cytometry), crucial control experiments should be performed to tease out the effects of the long incubation with the dye and the effects of UV radiation, with or without TNF. Do the dye and UV light perturb natural signaling dynamics or other physiology of the cells, such as cell cycle and apoptosis?

B) Importantly, the authors should compare the transcriptional response of cells with or without exposure to dye and UV light to verify that similar results are obtained under both conditions.

2) With regards to the quantification of nuclear GFP-RelA, the authors state that the nuclear to cytoplasmic intensity, termed NCI, is an "internally normalized measure". Although this quantity has been used previously to represent nuclear localization of NF-κB (Nelson DE et al. Science 2004 and subsequent publications from the same group), more groups now approximate the nuclear concentration of NF-κB with the background-adjusted mean nuclear intensity (Lee TK et al. 2009 Science Signaling; Sung MH et al. 2014 Science Signaling; Lee REC, et al. 2014 Molecular Cell; Kellogg RA and Tay S, 2015 Cell), while minimizing the photobleaching which decreases fluorescent signals. NCI is not fully internally normalized, because both the numerator and the denominator fluctuate over time due to biological effects, not only due to technical effects. Would the same damped oscillations dynamics be observed with background-adjusted mean nuclear intensity?

3) Regarding the comparison with the recent report by Kellogg and Tay (2015 Cell), one important difference that is not highlighted here is that the applied pulse characteristics are different. The Tay group delivered a pulse by an infusion of TNF-α without washing (i.e. T2 = 0, using the terminology of Figure 1A), assuming a clearance process through internalization and degradation. Zambrano et al. enforced a washing and reset period between consecutive pulses. While it is a matter of debate which method better mimics a physiologically relevant situation, the different forcing protocols may have affected their results and make a direct comparison problematic. Could this explain why the two groups observe different results?

4) Concerning the choice of pulse duration for different forcing frequencies, the authors should justify why different "T1" periods were used for different forcing frequencies (45 min for 90 min pulses, 30 min for 180 min pulses, etc.). Michael White's group has reported that the temporal profile of nuclear NF-κB depends on the pulse duration (Ashall L et al. 2009 Science, FigureS7). Fixing the pulse duration might simplify the interpretation of results from a systematic comparison.

5) With regards to the computational model, one concern is that many parameter values seem ill-defined experimentally. Although the model builds on a previously published version, a more exhaustive presentation of the model (to explain which parameters are fixed, and which are fitted) should be included. It is also recommended that experimental effort be made to narrow down the number of free parameters: this greatly improves the quality of biochemical models, and makes for better predictions. Finally, IκBα is encoded in the model to exert a "transcriptional repressor" role although the cited study does not present particularly convincing evidence for any direct repressor activity of IκBα. Is this term necessary to recapitulate the experimental data, especially for a minimal mathematical model?

6) With respect to the sensitivity of model behavior to parameter values, the authors acknowledge large parameter sensitivity of their model, with 2-fold changes in expression levels of proteins switching cells from damped to sustained oscillatory responses (c.f. Figure 1—figure supplement 5A-D).

A) An exact description of the parameters that were varied in these supplementary figures needs to be included.

B) More generally, a systematic parameter sensitivity of the model (e.g. in terms of generating varied types of oscillation) should be performed to allow readers to understand the robustness of the model (in terms of its predictive power). In particular, this could potentially identify key proteins whose up/down regulation would most affect the response phenotype. In turn, the natural variability in the expression levels of these signaling components should be taken into account: it would explain the phenotypic variability that must exist within the isogenic population of cells used in these experiments, but that was overlooked so far.

C) Most relevant to the biology of NFκB signaling, such parameter sensitivity would help reconcile the diverse dynamics of NFκB response depending on the cell type under consideration (c.f. comparison with Kellogg and Tay, 2015 study). For that reason, it would be particularly useful to document systematically the variety of phenotypes (from damped oscillations, to spontaneous oscillations etc.) to predict and explain how different cell types use NFκB response differently to generate different functional responses, and to support the authors' argument that damped oscillations responses may be common while the stable oscillations may only be a rare case.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "NF-κB oscillations translate into functionally related patterns of gene expression" for further consideration at eLife. Your revised article has been favorably evaluated by Aviv Regev (Senior editor) and a Guest Reviewing editor. The manuscript has been improved but there are some remaining minor issues that need to be addressed before acceptance, as outlined below.

In their revised submission, Zambrano and colleagues have addressed most of the concerns raised in the original review of this manuscript and the submission is accordingly much stronger. In particular, the authors have added a number of datasets to validate their experimental setup for live-cell imaging, mitigating initial concerns about damaging conditions for the cells under study. The new justification in the main text (Results) along with new supporting data (Figure 1 supplements) better explains the rationale for using this particular measure of nuclear NF-κB in this study. Furthermore, the detailed explanation now included in the Materials and methods presents a complete argumentation for this approach with their current datasetMaterials and methods, and clearly lays out the assumptions made and therefore provides an excellent overview for all readers familiar, or unfamiliar, with these types of data. Additional data and explanations in the text also better support the choice of stimulation profile. The new manuscript also provides a well-done analysis of their computational model which helps position their findings with regards to the dynamic behavior of TNF-induced NF-κB responses in a greater context of possible and probably behaviors. Only a few minor issues are noted for the revised version which we hope can be rapidly addressed by the authors.

Cell death, cell division and imaging conditions – overall, the additional data provided by the authors, particularly Figure 1—figure supplement 7 and Figure 5—figure supplement 2 do suggest that the cell death that is observed in some of the experiments is likely due to the presence of TNF (under flow condition and in the absence of serum), and not just an artefact of the use of a Hoescht dye and UV illumination to identify the nucleus of each cell. However, the evidence from PO4-gammaH2AX in Figure 1—figure supplement 6 has an important weakness and should probably be de-emphasized. Indeed, this epitope marks double strand breaks (DSBs) and while doxorubicin (positive control) is known to produce DNA double strand breaks (DSB) in cells in most phases of the cell cycle (DSB are linked to genes actively transcribed; e.g. http://www.ncbi.nlm.nih.gov/pubmed/25705119) UV illumination creates other types of DNA damage (pyrimidine dimers) but only lead to DSBs in cells in S-phase and on a longer time scale than that used bythe authors in this figure (6-12 hrs, vs. 3 hrs – e.g. http://www.ncbi.nlm.nih.gov/pubmed/20947453 and http://www.ncbi.nlm.nih.gov/pubmed/18800352). Therefore, lack of PO4-gammaH2AX staining does not necessarily preclude the possibility of UV-induced DNA damage.

Square vs. sawtooth stimulation profile. Because the geometry of the microfluidic chamber is very different in this study from that of the device used in Kellogg and Tay, it is unclear if a significant reduction in TNF concentration is achieved when medium flow is stopped (as noted by the authors, TNF is stable in their chamber, and would decay only if actively degraded by cells). Importantly, this could be an alternative explanation for the lack of synchrony – this stimulation profile may be closer to static concentration than to oscillatory stimulation – which should be noted in the discussion of these data.

Regarding the sensitivity analysis of the model – the sensitivity analysis does provide a greater context for the observations by the authors as well as the cited work by Kellogg and Tay. The authors mention that parameters related to the A20 feedback are the only one for which there is a clear difference between values that lead to different types of dynamics. When discussing this result, the authors may find it helpful to refer to Werner et al. (2008; Genes and Dev.), where the A20 feedback was explored computationally and experimentally. This paper should also be cited in the fourth paragraph of the Introduction, when introducing the A20 feedback.

Regarding the role of IκBα as a repressor – the new results from the authors do shed more light on that issue. However, the authors should review their model diagram as some are ambiguous or possibly misleading. It would be easier for the readers to interpret the model if the "activation" arrow for NF-κB and the "inhibition" arrow for IκBα were carefully positioned to target the reactions on which they act.

In the third paragraph of the subsection “Periodic forcing turns heterogeneous NF-κB oscillations into synchronous oscillations”:

a) "Replicate and show…"; "divide" or "proliferate" would be a better term to characterize the cells' behavior.

b) “Limited cell death… "; it may be more appropriate to just state "cell death…".

In the fourth paragraph of the Introduction: "In resting cells, p65:p50 exists mostly as a cytoplasmic complex bound to the IκB inhibitors (Hoffmann et al., 2002)" would be better rephrased to: "In resting cells, p65 exists mostly within a cytoplasmic complex bound to the IκB inhibitors (Hoffmann et al., 2002)."

"Figure 4: The authors should indicate in the main text or in the figure legend how many forcing cycles the cells were stimulated with before withdrawal of the forcing." The authors have now added cycle numbers on the figures although it would be helpful if the figure legend explained what these numbers represent.

https://doi.org/10.7554/eLife.09100.052

Author response

Summary:

How transcription factors affect responses to extracellular signals remains a key question in systems biology for which the NF-κB signaling pathway in mammalian cells has been a canonical study system. The manuscript "NF-κB oscillations translate into functionally related patterns" by Zambrano et al. presents a systematic investigation of single-cell NF-κB signaling dynamics and of population-level transcriptional responses in GFP-RelA knock-in mouse fibroblasts stimulated with controlled TNF pulses. The authors use a microfluidic device to deliver TNF to cells while avoiding the accumulation of secreted proteins. Although, this does not exclude cell-centric feedback regulation of the signaling response to TNF, the minimal ODE model proposed by Zambrano et al. nonetheless accounts well for their experimental observations and the authors conclude that NF-κB behaves as a damped oscillator in these fibroblasts. That said, the reviewers have raised serious concerns about the single-cell data (detailed below) that need to be addressed. The second half of the study, focused on the gene expression kinetics and functional correlates in cells treated with TNF pulses, leverages the use of oscillatory stimulations to reveal different clusters of genes that share dynamic patterns of expression. This is the most convincing and interesting finding of the manuscript. An important biological insight is that different functional pathways exhibit different expression dynamics; chemokine/chemokine receptors oscillate with TNF inputs, while immune responses integrate over longer timescales. The authors convincingly show that the decay rate of the mRNA of a gene of interest is sufficient to account for the different dynamic regimes of gene regulation. We thank the reviewers for the insightful comments, which highlight the need for important controls that were lacking or had not been mentioned in the previous version. We describe below how we have addressed the major concerns of the reviewers about the single-cell data and we describe the changes introduced in the manuscript.

Our controls reinforce the conclusion that NF–κB oscillations can lock to a wide variety of input stimuli as damped oscillations. We further argue that our isogenic cell population behave as heterogeneous damped oscillators for a constant stimulus, and this is supported by further imaging experiments and mathematical modelling.

Essential revisions: 1) The live cell imaging conditions seem to have been very damaging to cell physiology. Stress from photo-damage during live-cell imaging is a very serious issue which has hampered many cell biological studies (Cole R, 2014 Cell Adhesion & Migration 8 (5) 452-259), and it is essential to address this problem. In the data presented by Zambrano et al., the supplementary movies show that the knock-in fibroblasts express very low levels of GFP-RelA making imaging challenging and it is quite evident from these movies that there is significant cell death and stress, and few cell divisions. Photo-stress was likely increased in this study due to the repeated exposure to the UV laser used for Hoechst imaging. A key concern is that with such a high level of stress, the observed NF-κB dynamics may not represent the true signaling dynamics induced solely by TNF. This could explain the authors observe highly damped oscillations, which are different from the sustained oscillations observed and reported in their first paper on the same GFP-RelA knock-in MEFs (Sung MH et al. 2009 PLoS One). In contrast to the manuscript under consideration, no UV or Hoechst were used in the 2009 study and healthy cell divisions were frequently observed throughout the long-term imaging supplementary movies. In light of these important concerns, the inclusion of the single-cell imaging data in the paper requires additional control data: A) Because nuclear dyes can alter cell cycle and/or cause oxidative stress and mutations (Martin RM, Leonhardt H, and Cardoso MC, 2005 Cytometry; Ge J, Wood DK, et al. 2013 Cytometry), crucial control experiments should be performed to tease out the effects of the long incubation with the dye and the effects of UV radiation, with or without TNF. Do the dye and UV light perturb natural signaling dynamics or other physiology of the cells, such as cell cycle and apoptosis?

On one hand, these comments highlight the need to complete and make explicit the controls already made that exclude any interference of our nuclear labeling with the natural signaling of NF-κB in the cells. On the other hand, they suggest that our description of single cell dynamics as damped (and not sustained) oscillations was not as sharp as it should have been.

To address the first issue, we first performed several control experiments.

DNA damage:

To exclude DNA damage in the cells due to imaging, we assessed by immunostaining the presence of H2AX phosphorylation (gammaH2AX), a marker for ongoing DNA damage repair. We compared gammaH2AX staining in cells that were not exposed to Hoechst nor imaged and in those exposed to Hoechst and/or laser irradiation for imaging.

Figure 1—figure supplement 6 shows the results of this analysis: no difference in gammaH2AX intensity was detected in the cell nuclei upon different treatments, only the expected basal levels are observed (Turinetto & Giachino, 2015). It is also worth considering that NF-κB plays a pivotal role in activating the DNA damage response upon irradiation and a nuclear translocation would be expected if imaging is anyhow harmful. Importantly, we do not see increasing levels of global activation of NF-κB through the NCI time series for unstimulated cells (Figure 1—figure supplement 8A and Video 6).

Hence, at the low concentration of labeling used and the relatively low energy of our laser (405 nm) DNA damage remains almost undetectable. These results have been discussed in the main text and described in detail in the Materials and methods section.

Cell divisions:

It is important to point out that cells are cultured in medium containing reduced serum (0.1% FBS) to minimize background noise during imaging, a necessary technical condition when the signal to noise ratio is very low. In low serum, cells progressively transit through the G2/M phase and accumulate in G1, as documented by the cell cycle analyses reported in Figure 1—figure supplement 10. The low serum does not block the cell cycle completely and cells in mitosis are present until the very last imaging frame in all conditions tested. Cells are healthy as can be appreciated in the movies of unstimulated cells (Video 6). This is the main explanation for the low number of mitoses noted by the reviewers, whom we thank for asking for clarifications.

We acknowledge that a more accurate assessment of the imaging effect on the cell cycle would have been provided by single-cell based approaches like the FUCCI cell cycle sensor (https://www.thermofisher.com/it/en/home/life-science/cell-analysis/cell-viability-and-regulation/cell-cycle/live-cell-imaging-of-cell-cycle-and-division.html) that allows following cell divisions within a cell population during imaging. We actually tried this approach but our fibroblasts are very resistant to transfection. The minimal fraction of cells that were labeled with the sensor did not allow a statistically significant evaluation of cell cycle alteration, if any, in the cells undergoing imaging in the presence/absence of TNF-α.

Cell death:

In our previous work with these cells and imaging settings (Zambrano et al., 2014a) we quantified photo-stress after 15 hours of imaging of Hoechst-exposed (Video S1 and S2 in the paper) and non-exposed cells (Videos S3 and S4 of that paper), finding equally low death rates (close to 5%). Indeed, in microfluidics set-up, cell death is even lower in unstimulated cells than in Zambrano et al. 2014, see Video 6.

However, the reviewers are right in pointing out that Video 1 with cells treated with 10ng/ml TNF-α shows considerable cell death and now we propose an explanation for that effect. We suspected that TNF-α–but not Hoechst/UV imaging– might lead to cell death. To test this we evaluated the appearance of apoptotic cells in cell populations exposed to a continuous flow of different doses of TNF-α in the microfluidics chamber without Hoechst/UV imaging. In Figure 1—figure supplement 7 we show GFP-imaging fields from time-lapse acquisitions for 0.1, 1, and 10 ng/ml TNF-α continuous stimulation. The appearance of dead cells is considerably anticipated as the dose of TNF-α increases: the first dead cells appear on average at 12, 9 and 3 hours, respectively (red frames).

FACS analysis of the cell cycle supports this conclusion while excluding Hoechst toxicity at population level. Untreated cells or cells exposed to Hoechst, in the presence or absence of 10ng/ml of TNF, showed no global perturbation of the cell cycle but apoptotic cells appeared after three hours of incubation with TNF-α (Figure 1—figure supplement 10). Interestingly, a similar phenomenon was also observed by the Tay group: cells under static stimulation with 10 ng/ml TNF-α (movie S1 in (Tay et al., 2010) undergo almost no cell death, while many cells under constant TNF-α flow die (movie S1 in new Kellogg & Tay, 2015). This probably depends on culture conditions and on the cell type; the precise mechanism will be the object of future investigation.

Sustained versus damped oscillations:

The reviewers are concerned that in our culture conditions oscillations were “damped” and thus different from the “sustained” oscillations that were observed for the same cells in Sung et al. 2009 using different culture conditions. This is in part a semantic problem involving the terminology used. In Sung et al. 2009 oscillations of GFP-p65 were referred to as “sustained” as opposed to the extremely damped dynamics of NF-κB revealed by biochemical population-level assays as in the seminal paper of Hoffmann et al., 2002.

In the present work we do not claim that our cells have different dynamics from the one described in Sung et al. 2009 and Zambrano et al. 2014. We described the dynamics as “damped” as opposed to “sustained oscillations”, the latter referring to the regular series of peaks – almost sinusoidal – for time lapses longer than 12 hours reported by Kellogg and Tay, 2015.

In dynamical systems terms, damped oscillations are those that tend to stabilize to a fixed point while sustained oscillations converge to limit cycle. This is what we mean by damped oscillations here. In Figure 1—figure supplement 2D, we show a representative collection of the heterogeneous dynamics observed under continuous flow of TNF-α, from non-oscillating cells to increasingly oscillating cells, as previously reported. This heterogeneity was summarized in the number of peaks per cell shown in Figure 1—figure supplement 3B. This is remarkable provided that culture conditions are not static as in (Sung et al., 2009; Zambrano et al., 2014a), so autocrine-paracrine signaling, which might induce oscillations in late phases of the stimulation, is not possible.

Overall, though, the dynamics we observe cannot be described as sustained using Kellog and Tay terminology, provided that the peaks we find become more irregular and infrequent as the stimulation proceeds, in agreement with what observed previously. Most of our time series oscillate to apparently end up fluctuating around a value higher than the pre-stimulation level. In conclusion, we unfortunately used the same term “sustained” in two different meanings in two different papers, but there is little we can do: we should retroactively change the term used in Sung et al. 2009. We now try to define exactly in our manuscript the term sustained as we use it.

Possible interference by Hoechst and imaging on oscillations:

Finally, to further assess if labeling and imaging induce distortions of NF-κB dynamics, we decided to manually quantify NF-κB dynamics for cells that were not exposed to Hoechst nor imaged with UV. Examples of the manually derived single cell traces are shown in Figure 1—figure supplement 8C, which show a heterogeneous but damped oscillatory dynamics as the ones shown in this manuscript and in our previous works, under imaging in the presence of Hoechst.

Actions taken:

In order to sharpen the manuscript, we have briefly delineated in the section “Periodic forcing turns heterogeneous NF-κB oscillations into synchronous oscillations” that controls were performed to exclude any interference between the signalling dynamics of NF-κB and the imaging conditions, and that cell death is mostly due to the continuous flow of TNF-α, referring to the movies. This is illustrated with several supplemental figures cited in the main text and thoroughly described in “Safety assessment for Hoechst staining” in the Materials and methods section.

To clarify further our description of the oscillations as damped, in section “Periodic forcing turns damped NF-κB oscillations into synchronous oscillations” we refer to the new figure, Figure 1—figure supplement 2D, to illustrate the heterogeneous damped oscillatory dynamics that we observe as opposed to the sustained oscillations that are observed by Tay’s group. We also put our results in the context of previous observations done with the same cells.

B) Importantly, the authors should compare the transcriptional response of cells with or without exposure to dye and UV light to verify that similar results are obtained under both conditions.

This is actually an important point that needed clarification in the manuscript. In fact, we did perform our experiments of transcription with cells that were not exposed to Hoechst nor UV illuminated. As shown by the controls performed to address comment 1A, it is clear that the dynamics are not distorted by the imaging conditions, so neither should be the transcriptional response. However, in order to further exclude effects on transcription, we performed RT-PCR transcription assays for IκBα and Ccl5 in cells under constant TNF-α both with and without Hoechst staining. The results shown in Figure 5—figure supplement 2 indicate no significant difference between the two conditions.

Actions taken:

This is now clarified in the revised version of the manuscript. We also added in Materials and methods that no nuclear labelling was used in the transcription experiments, and that we performed controls showing that nuclear labelling wouldn’t have altered the results anyway. There we refer to the newly added Figure 5—figure supplement 2.

2) With regards to the quantification of nuclear GFP-RelA, the authors state that the nuclear to cytoplasmic intensity, termed NCI, is an "internally normalized measure". Although this quantity has been used previously to represent nuclear localization of NF-κB (Nelson DE et al. Science 2004 and subsequent publications from the same group), more groups now approximate the nuclear concentration of NF-κB with the background-adjusted mean nuclear intensity (Lee TK et al.

2009 Science Signaling; Sung MH et al. 2014 Science Signaling; Lee REC, et al. 2014 Molecular Cell; Kellogg RA and Tay S, 2015 Cell), while minimizing the photobleaching which decreases fluorescent signals. NCI is not fully internally normalized, because both the numerator and the denominator fluctuate over time due to biological effects, not only due to technical effects. Would the same damped oscillations dynamics be observed with background-adjusted mean nuclear intensity?

This comment made us realize that we did not provide convincing enough arguments of the validity of Nuclear to Cytoplasmic Intensity (NCI) as a quantifier of NF-κB dynamics and on the rationale of our selection. To summarize, we believe that the question raised by the reviewers can be answered with a clear yes and detail our reasoning below.

First, the imaging rationale: As we can see in some of the movies accompanying this manuscript (see e.g. Videos 1 and 6) there are global fluctuations in the intensity of the images and of the background, which also vary in space. This is essentially due to small changes in the focus, slight inclination of the plates and changes in the laser intensity that might amplify-or reduce- the amount of light measured. These changes are sizeable when the signal to noise ratio is very low, as in the case of physiological levels of p65. Photobleaching, if any, is negligible: one can see in our movies – as in Video 1 – that cells are bright even in late phases of the imaging. Thus our rationale is to use a ratio of background-corrected intensities, so that these distortions will cancel out; this is the same rationale that we applied in our previous works, where we used ratios of intensities to quantify the dynamics (Sung et al., 2009; Zambrano et al., 2014a).

Second, the mathematical rationale: A precise description of our rationale for selecting NCI has been reported in the extended paragraph “Quantification of NF-κB dynamics” in the Materials and methods section. The mathematical discussion leads to conclude that, being a ratio of intensities, NCI is a faithful indicator of the ratio between nuclear and cytoplasmic NF-κB and hence robust to global changes in the image intensity. We use the same reasoning to show why the background-corrected average nuclear intensity would reflect the global imaging distortions in our settings and hence would not be an adequate readout of the nuclear amount of NF-κB. However, as the reviewers rightly point out, the numerator and denominator of this ratio can fluctuate due to biological reasons and this might blur or overestimate the oscillatory behavior of nuclear NF-κB, which is the truly important readout. If it were possible to assess that total NF-κB is constant, then the oscillations in the nuclear NF-κB would correspond to oscillations in NCI, as we demonstrate in Supplemental Information, due to the fact that the latter is a monotonously increasing function of the former, (Figure 1—figure supplement 8B). We then provide evidence that NF-κB is indeed constant (see below).

Third, the biochemical rationale: Is the total amount of NF-κBa constant parameter? Although this has been widely assumed in the field and included in different mathematical models present in the literature, from the seminal paper (Hoffmann et al., 2002) to more recent papers as (Kellogg & Tay, 2015), this view has been challenged by a recent work (Sung et al., 2014) showing that p65 is regulated by a positive feedback in macrophages under LPS stimulation. Thus, to confirm that the “constant p65” assumption holds true for our GFP-p65 MEFs, we quantified p65 in cells stimulated with 10 ng/ml TNF-α for 8 hours with sampling every hour. Immunoblots in Figure 1—figure supplement 9 show no large changes in the p65 levels upon stimulation. We did not quantify p65 at single cell level using microscopy because fluctuations of our experimental setup might lead to imprecise results. However, we consider our assumption valid and in line with the literature.

Finally, variations in the areas of the nucleus and the cytoplasm might introduce small distortions, although they typically vary slowly. This problem can be easily solved by discarding cells for which the areas change abruptly – an indicator of imminent mitosis or cell death –, as our software does.

To further support the validity of NCI we are providing additional data.

First, we have that NCI(t) remains constant at single cell level – except for some spontaneous oscillations, as previously reported (Zambrano et al., 2014a) – as shown in Figure 1—figure supplement 8A. This is remarkable because the movie from which these series were obtained shows considerable decrease in image intensity (see Video 6). Along the same lines we do not find upwards or downwards average trends in our data of NCI time series for different conditions, which indicates further that the time series are properly normalized and reflect well the dynamics. We have also plotted side-by-side the NCI values and average nuclear intensity obtained by manual segmentation in Figure 1—figure supplement 8C. For this experiment the imaging conditions were stable so the same kind of qualitative behavior is obtained, as expected. A last indicator of the sensitivity of the quantifier comes from its ability to detect even small oscillatory peaks, see Figure 3 and Figure 3—figure supplement 1 as an example for 0.1 ng/ml TNF-α.

Actions taken:

We have tried to make all this clearer in the revised manuscript. In the first section of the results, where we briefly describe our approach, we explain that NCI is a valid quantifier of NF-κB dynamics while being robust to imaging distortions. Also, we decided to abandon the “internally normalized” term, which might be misleading.

Finally, in the Materials and methods section, we have introduced a new paragraph entitled “Selection of the quantifier of NF-κB dynamics” where we reproduced the discussion above on the validity of each quantifier providing the appropriate mathematical demonstration. The cited movie and supplement figures have been included in the new manuscript.

3) Regarding the comparison with the recent report by Kellogg and Tay (2015 Cell), one important difference that is not highlighted here is that the applied pulse characteristics are different. The Tay group delivered a pulse by an infusion of TNF-α

without washing (i.e. T2 = 0, using the terminology of Figure 1A), assuming a clearance process through internalization and degradation. Zambrano et al. enforced a washing and reset period between consecutive pulses. While it is a matter of debate which method better mimics a physiologically relevant situation, the different forcing protocols may have affected their results and make a direct comparison problematic. Could this explain why the two groups observe different results?

We thank the reviewers for giving us the opportunity to expand on this point, which now has been addressed in the revised manuscript, also by adding new experimental results.

Indeed a crucial difference of the setup of Kellogg and Tay with respect to ours is that they periodically and quickly refresh the medium in contact with the cells, after which they assume that a clearance process takes place. This should produce a sawtooth TNF-α profile. We, instead, do constantly control the amount of TNF-α in direct contact in the cell by subjecting them to a continuous flow of the selected doses D1 and D2 for different times T1 and T2. Importantly, as we already pointed out, our approach also should reduce to the minimum the contribution of autocrine-paracrine signaling in the dynamics.

To gain further insights in whether these differences could explain the observed differences in the synchronization mechanism, we have measured the dynamics of our cells under sawtooth-like profiles of TNF-α with different periodicities. The results are shown in the new Figure 4—figure supplement 4. GFP-p65 MEFs synchronize to sawtooth pulses with Tf=90'. In contrast to Kellogs and Tay’s results, our MEFs synchronize also to Tf=60'. However, the synchronization after the first pulse is quickly lost for Tf=180'.

Interestingly, with the sawtooth approach synchronizations are less sharp and reminiscent of those obtained using alternating doses (D1>0 and D2>0) for Tf=60' (Figure 4—figure supplement 1B-D), Tf=90' (Figure 2—figure supplement 3) and Tf=180' (Figure 3—figure supplement 2). These observations suggest that culturing conditions, absence of flow and chamber geometries possibly affect TNF-α decay and clearance. Moreover, autocrine-paracrine signaling may further contribute to distortions of the dynamics.

However, coming back to the point raised by the reviewers, the shown ability of our cells to synchronize to a sawtooth of Tf=60', contrarily to what is reported by Kellogg and Tay (2015 Cell) for a similar forcing, further suggests that the synchronization mechanism strictly depends on cellular and environmental factors.

Actions taken:

We have included the new Figure 4—figure supplement 4 in the manuscript and discussed it in the section The synchronization mechanism of GFP-p65 knock-in MEFs is not entrainment. We have put this figure in the context of our discussion, stating that it provides further confirmation to our interpretation of the synchronization mechanism as that of a damped oscillator to a (sufficiently) strong external periodic forcing.

4) Concerning the choice of pulse duration for different forcing frequencies, the authors should justify why different "T1" periods were used for different forcing frequencies (45 min for 90 min pulses, 30 min for 180 min pulses, etc.). Michael White's group has reported that the temporal profile of nuclear NF-κB depends on the pulse duration (Ashall L et al. 2009 Science, FigureS7). Fixing the pulse duration might simplify the interpretation of results from a systematic comparison.

We agree with the reviewers: using the same value of in all the experiments would help in the interpretation of the results. A comprehensive analysis of many different combinations of T1, T2, D1 and D2 would have been extremely time-consuming and economically challenging. Therefore we started exploring only few combinations. We discuss here the reasons of our selection of T1 in each experiment while providing additional supplementary experiments to facilitate data comparison, and reinforcing the central message of the dynamics part of our work: that NF–κB oscillations can lock to a wide variety of input stimuli. This message has been highlighted further in the manuscript.

In our first experiments with Tf =90 min we decided to use T1 =45 min to see whether a "long square wave" forcing could give rise to synchronous dynamics similar to the T1 =5 min short pulses documented in the paper by Michael White's group (Ashall et al., 2009). To explore lower frequencies, like Tf =45 min, smaller values of T1have been chosen. We first tried with T1 =30 min, T2=15 min, (Figure 3—figure supplement3, upper panels) and found that synchrony to the forcing was very poor, which highlighted the importance of a resetting phase in order to obtain synchronization to the input. By reducing the value of T1 =30 min to the value of T1 =22,5 min shown in Figure 3, we find the synchrony of the oscillations for different doses.

On the other hand, we selected T1=30 min for the experiment with Tf =180 min shown in Figure 3 because we observed that T1=30 min for Tf =90 min produced synchronous dynamics (Figure 3—figure supplement 3, lower panels), almost as clear as the ones observed for T1 =45 min. We also speculated that T1 =30 min,Tf =180 min might produce even sharper oscillations in the transcriptional output, which in turn was useful for Figures 5 and 6.

Actions taken:

We have better motivated our selection of T1in the manuscript and have included Figure 3—figure supplement 3 where the interested reader can find graphs on the dynamics for the conditions described above. Furthermore, we have highlighted the importance of the resetting time in obtaining the synchrony that we observe in the Tf =45 min, which is something that we overlooked in the original version of our manuscript.

5) With regards to the computational model, one concern is that many parameter values seem ill-defined experimentally. Although the model builds on a previously published version, a more exhaustive presentation of the model (to explain which parameters are fixed, and which are fitted) should be included. It is also recommended that experimental effort be made to narrow down the number of free parameters: this greatly improves the quality of biochemical models, and makes for better predictions. Finally, IκBα is encoded in the model to exert a "transcriptional repressor" role although the cited study does not present particularly convincing evidence for any direct repressor activity of IκBα. Is this term necessary to recapitulate the experimental data, especially for a minimal mathematical model?

Model description:

We agree with the reviewers that many parameters of the model are ill-defined: there is always a trade-off between the complexity of the model designed and the number of parameters that will be well defined, provided that any simplification implies to conflate in few processes many different ones. This is the case of our 8-dimensional mathematical model, but this holds true also for higher dimensional mathematical models of the NF-κB system, as e.g. Tay et al., 2010 and Ashall et al., 2009 – with twice as many variables as ours – where still ill-defined parameters can be found. As the reviewers implicitly suggest, a clear description of how the model is built can help to find a correspondence between variations in the experimental conditions and the system’s parameters. Hence, we have introduced in the Materials and methods a detailed description of the processes that our mathematical model takes into consideration, reproducing the basic discussion of the work where we described a previous version of the model (Zambrano et al. 2014b), including the normalizations performed, which leads to a simplification and a reduction of the free parameters.

Parameter variation:

Along the same lines, we have followed what we did in (Zambrano et al., 2014b) and assigned to each of the model parameters an uncertainty degree D, which is low D=0.2 for parameters derived from biochemical rates of the literature (in our case, mostly from Tay et al., 2010, such as decay rates of proteins and RNA) and high (D=1) for others that are unknown or that correspond to processes that have been simplified, as in our case those involved in the A20 negative feedback. Parameters are then allowed to be varied by a random factor in the [10-D, 10D] interval, which means that the largest uncertainty degree imply a variation of up to one order of magnitude around the initial parameter value. In the new Table 2 we now explicitly state the biochemical rates and constants used, their origin, the parameters of the model to which they are related and the uncertainty degree considered for each of them. Importantly, in our fittings we allow each of the parameters to vary only by a factor in the [10-D, 10D] interval. This is now clearly stated in the Materials and methods.

Actions taken:

We also have clarified the notation of the parameters making explicit those that are fixed and those that are fitted. We essentially have three sets of parameters: PS, PNF-κB, and PG. The first, PS, accounts for the parameters of the external time-varying signal; PNF-κB, is the set of parameters that regulate NF-κB dynamics via the double IκBα–A20 negative feedback loop, and PG encloses the parameters used to reproduce the gene expression. Thus, the mathematical description of the expression of two genes under the control of NF-κB under the same external forcing conditions would have the same values of PS and PNF-κB, but different –individually fitted – values of PG. In Figure 5, then, the same values of PNF-κB where used to reproduce the NCI, IκBα and CCL5 profiles, but for each of these a different PG was used. Similarly, in the fittings of Figure 6 we use the same values of PS and PNF-κBas in Figure 5 but then each single gene is fitted with its own PG. This is also clearly explained in the new version of the manuscript, as well as in the captions of Figures 5 and 6.

On the role of IκBα as a repressor:

Finally, concerning the role of IκBα as a repressor, we agree that the evidence provided in (Arenzana-Seisdedos et al., 1995) is not conclusive. However, the role of IκBα as a repressor has been considered in a number of works, e.g. (Kellogg & Tay, 2015; Lipniacki et al., 2006; Tay et al., 2010), while not in others (Ashall et al., 2009; Sung et al., 2014). We have tested if our model can reproduce observed profiles even if such transcriptional repressor is not present, when simulating both the dynamics of NF-κB (Figure 6—figure supplement 7A) and gene expression (Figure 6—figure supplement 7B). Gene expression profiles are indeed reproduced, for genes with both increased and decreased transcription, with relative errors comparable to those obtained with our original model (Figure 6—figure supplement 7D). It can also be observed that the degradation rate is the key parameter to reproduce the dynamics of each cluster of genes (Figure 6—figure supplement 7E). Importantly, this exploration made us realize that for clusters 4-6 containing genes with decreased transcription we were essentially fitting the gene expression profiles as a simple RNA degradation process.

Actions taken:

All this has now been clarified in the text, where we also mention the newly added Figure 6—figure supplement 7. We also clarify that the role of IκBα as a repressor is not fundamental.

6) With respect to the sensitivity of model behavior to parameter values, the authors acknowledge large parameter sensitivity of their model, with 2-fold changes in expression levels of proteins switching cells from damped to sustained oscillatory responses (c.f. Figure 1—figure supplement 5A-D). A) An exact description of the parameters that were varied in these supplementary figures needs to be included.

Using the notation that we have introduced to address point number 5, we have now specified that the parameters varied were PNF-κB, according to the uncertainty degree that we assigned to each of them. We have added a similarly precise description of the parameters varied in the remaining numerical supplemental figures.

B) More generally, a systematic parameter sensitivity of the model (e.g. in terms of generating varied types of oscillation) should be performed to allow readers to understand the robustness of the model (in terms of its predictive power). In particular, this could potentially identify key proteins whose up/down regulation would most affect the response phenotype. In turn, the natural variability in the expression levels of these signaling components should be taken into account: it would explain the phenotypic variability that must exist within the isogenic population of cells used in these experiments, but that was overlooked so far.

Following the reviewers’ advice, we have performed a numerical exploration of our model in order to identify the variety of oscillatory behaviours that it can display. We varied the parameters of the model according to the uncertainty degree assigned to them. Since our model is a system of ordinary differential equations, it is possible to determine numerically the stability of the fixed points of the flow determined by equations (1)-(8) (see Materials and methods), and thus what dynamics should arise for each parameter combination.

We did generate a library of 10,000 random parameter combinations in the intervals prescribed, and did pre-select those that gave a response in the time series of N(t). The eigenvalues of the fixed point are plotted in the complex plane in Figure 1—figure supplement 5A; those in red correspond to sets of eigenvalues for which at least one of them has a real positive part. In that situation, the system converges to sustained oscillations, but the fraction of those combinations is less than 10% of the total, as shown in Figure 1—figure supplement 5B. Examples of sustained oscillations similar to the ones shown by Kellogg and Tay 2015 and some spiky oscillation, similar to the ones shown in Ashall 2009, are shown in Figure 1—figure supplement 5D.

We also find that most of the trajectories have four eigenvalues with imaginary parts, see Figure 1—figure supplement 5C; this implies that most of the trajectories have two competing frequencies. For this reason, we find examples in which the interaction of fast and slow oscillating timescales can be observed. Additional examples of the variety of time series present can be shown in panel Figure 1—figure supplement 5D. The fact that the signaling of NF-κB can present rich dynamics was pointed out in Sung et al. 2009. We hope that our results will contribute to further illustrate this important point that shows how even in isogenic populations heterogeneous dynamics – due to different parameter combinations – might arise.

Finally we present in Figure 1—figure supplement 5C the intervals of the parameters that give rise to sustained or damped oscillations (recall that the parameters selected are those giving trajectories responsive to the stimulus). We find that in our simple model the parameter ranges are very similar in oscillating and non-oscillating trajectories. Only for certain parameters related to the A20 feedback (n and dA) we find a clear difference. This suggests that variations in this upstream negative feedback can lead to strong variations in the type of dynamics. Overall, though, our results suggest that only the precise combination of parameters is what determines whether the system displays sustained or damped oscillations. Only further analytical work on the structure of bifurcations in parameter space may shed light on this important point. This type of work could also pave the way for targeted experiments that might be able to identify key proteins in NF-κB regulation, but this is far from the scope of the present work.

Actions taken:

We have tried to make all this clearer in the revised manuscript. An appropriate description has been added to the first section of the Results. The details of the stability analysis have been added to the Materials and methods section. Figure 1—figure supplement 5 has been improved with panels summarizing the discussion above.

C) Most relevant to the biology of NFκB signaling, such parameter sensitivity would help reconcile the diverse dynamics of NFκB response depending on the cell type under consideration (c.f. comparison with Kellogg and Tay, 2015 study). For that reason, it would be particularly useful to document systematically the variety of phenotypes (from damped oscillations, to spontaneous oscillations etc.) to predict and explain how different cell types use NFκB response differently to generate different functional responses, and to support the authors' argument that damped oscillations responses may be common while the stable oscillations may only be a rare case.

We believe that some the results previously shown and the revised version of the manuscript can address this point. In response to Point 2 we generated the new Figure 1—figure supplement 2, where we document the variety of phenotypes that we observe with examples of trajectories with different degrees of dampening, although the majority of them are those with few peaks (as shown in Figure 1—figure supplement 3). This indeed illustrates that in an isogenic population different oscillatory phenotypes can be observed. Hence, as the reviewers point out, it should be expected that also different cell types will display different oscillatory behaviors: this is probably the case of those studied in the work by Kellogg and Tay 2015, but the same could applied to recently published results with fibroblasts (see e.g. Figure 5, (Cheng et al., 2015)) or macrophages (see e.g. Figure 4B in Sung et al., 2014), that present also remarkable variability of behaviors. Furthermore, the numerical exploration performed to address the previous point illustrates how variations in the parameter combinations can give rise to widely different dynamics, which hints in the same direction. We can only claim that in our cells damped oscillations are much more common; our model suggests that damped oscillations are more frequent. Only the study of a wide variety of cell types can tell which of the two types of dynamics is more frequent. This will be hopefully the direction of our future work.

Actions taken:

All this has been made clearer in the first Results section and in the Discussion, where we state that our results (numerical and experimental) can also provide clues to reconcile the visions emerging from works where different cell types were considered, that we cite in the new version.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

[…] Only a few minor issues are noted for the revised version which we hope can be rapidly addressed by the authors. Cell death, cell division and imaging conditionsoverall, the additional data provided by the authors, particularly Figure 1—figure supplement 7 and Figure 5—figure supplement 2 do suggest that the cell death that is observed in some of the experiments is likely due to the presence of TNF (under flow condition and in the absence of serum), and not just an artefact of the use of a Hoescht dye and UV illumination to identify the nucleus of each cell. However, the evidence from PO4-gammaH2AX in Figure 1—figure supplement 6 has an important weakness and should probably be de-emphasized. Indeed, this epitope marks double strand breaks (DSBs) and while doxorubicin (positive control) is known to produce DNA double strand breaks (DSB) in cells in most phases of the cell cycle (DSB are linked to genes actively transcribed; e.g. http://www.ncbi.nlm.nih.gov/pubmed/25705119) UV illumination creates other types of DNA damage (pyrimidine dimers) but only lead to DSBs in cells in S-phase and on a longer time scale than that used bythe authors in this figure (6-12 hrs, vs. 3 hrse.g. http://www.ncbi.nlm.nih.gov/pubmed/20947453 and http://www.ncbi.nlm.nih.gov/pubmed/18800352). Therefore, lack of PO4-gammaH2AX staining does not necessarily preclude the possibility of UV-induced DNA damage.

To address this, we assessed the presence of thymine dimers as a specific indicator of DNA damage upon UV exposure. Thymine dimers in imaged cells were at a level below the sensitivity threshold of the immunostaining method. Now results are reported in a new figure (Figure 1—figure supplement 6).

We agree with the reviewers’ comment on the role of gammaH2AX. The text in the Materials and methods has been changed to de-emphasise the relevance of this marker for photo-damage. Appropriate new references have been added.

Square vs. sawtooth stimulation profile. Because the geometry of the microfluidic chamber is very different in this study from that of the device used in Kellogg and Tay, it is unclear if a significant reduction in TNF concentration is achieved when medium flow is stopped (as noted by the authors, TNF is stable in their chamber, and would decay only if actively degraded by cells). Importantly, this could be an alternative explanation for the lack of synchronythis stimulation profile may be closer to static concentration than to oscillatory stimulation

which should be noted in the discussion of these data.

This is an insightful comment. Indeed, the different geometry of the chamber and presumably different culture conditions can lead to a different decay rates for the TNF-α when put in contact with the cells. This observation has been included in the discussion of the sawtooth profile data.

Regarding the sensitivity analysis of the model

the sensitivity analysis does provide a greater context for the observations by the authors as well as the cited work by Kellogg and Tay. The authors mention that parameters related to the A20 feedback are the only one for which there is a clear difference between values that lead to different types of dynamics. When discussing this result, the authors may find it helpful to refer to Werner et al. (2008; Genes and Dev.), where the A20 feedback was explored computationally and experimentally. This paper should also be cited in the fourth paragraph of the Introduction, when introducing the A20 feedback.

We thank the reviewers for this comment. Indeed, the paper mentioned is one of the first descriptions of how the interplay between the A20 and IκB regulatory modules provides a fine-tuning of NF-κB response to an external stimulus. We now cite it when briefly describing the A20 feedback in the introduction. We also cite it when discussing the role of A20 in the type of dynamics observed in the paragraph “Classification of the different dynamical responses for constant stimulus”of the“Mathematical modelling” section of the Materials and methods. Finally, we also cite it in the Discussion when commenting on the same results.

Regarding the role of IκBα as a repressor

the new results from the authors do shed more light on that issue. However, the authors should review their model diagram as some are ambiguous or possibly misleading. It would be easier for the readers to interpret the model if the "activation" arrow for NF-κB and the "inhibition" arrow for IκBα were carefully positioned to target the reactions on which they act.

The diagrams in Figure 5, Figure 1—figure supplement 4 and Figure 6—figure supplement 7 have been changed as suggested. Furthermore, the figure captions of these figures now state clearly the meaning of the green and red arrows (transcriptional activation and repression, respectively).

In the third paragraph of the subsection “Periodic forcing turns heterogeneous NF-κB oscillations into synchronous oscillations”:

a) "Replicate and show…"; "divide" or "proliferate" would be a better term to characterize the cells' behavior.

We have changed the text accordingly.

b) “Limited cell death… "; it may be more appropriate to just state "cell death…".

The text has been changed.

In the fourth paragraph of the Introduction: "In resting cells, p65:p50 exists mostly as a cytoplasmic complex bound to the IκB inhibitors (Hoffmann et al., 2002)" would be better rephrased to: "In resting cells, p65 exists mostly within a cytoplasmic complex bound to the IκB inhibitors (Hoffmann et al., 2002)."

This has been rephrased as suggested.

"Figure 4: The authors should indicate in the main text or in the figure legend how many forcing cycles the cells were stimulated with before withdrawal of the forcing." The authors have now added cycle numbers on the figures although it would be helpful if the figure legend explained what these numbers represent.

Now we describe in the legend of Figure 4 what those numbers represent.

https://doi.org/10.7554/eLife.09100.053

Article and author information

Author details

  1. Samuel Zambrano

    1. Division of Genetics and Cell Biology, San Raffaele Scientific Institute, Milan, Italy
    2. San Raffaele University, Milan, Italy
    Contribution
    SZ, Performed the experiments, Mathematical modelling and model fitting, Analysis of the dynamics, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  2. Ilario De Toma

    San Raffaele University, Milan, Italy
    Contribution
    IDT, Performed the experiments, Analysis of the dynamics, Bioinformatics analysis, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
  3. Arianna Piffer

    San Raffaele University, Milan, Italy
    Contribution
    AP, Performed the experiments, Mathematical modelling and model fitting, Analysis of the dynamics, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
  4. Marco E Bianchi

    1. Division of Genetics and Cell Biology, San Raffaele Scientific Institute, Milan, Italy
    2. San Raffaele University, Milan, Italy
    Contribution
    MEB, Analysis of the dynamics, Conception and design, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    bianchi.marco@hsr.it
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0002-5329-6445
  5. Alessandra Agresti

    Division of Genetics and Cell Biology, San Raffaele Scientific Institute, Milan, Italy
    Contribution
    AA, Performed the experiments, Analysis of the dynamics, Bioinformatics analysis, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    agresti.alessandra@hsr.it
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon 0000-0001-5006-9506

Funding

European Marie Curie Action 2011 (298447-NonLinKB)

  • Samuel Zambrano

Intra European Fellowships (MCA-2011-298447NonLinKB)

  • Samuel Zambrano

AIRC2010/2013 project R0444

  • Marco E Bianchi

Flagship program EPIGEN

  • Marco E Bianchi

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

SZ was supported by the Intra-European Fellowships for career development-2011–298447NonLinKB. The work was supported by the IEF-2011-298447 research funding to SZ, the Italian AIRC 2010/2013 project R0444 and Flagship Program EPIGEN to MEB. We are grateful to Rohini Nair and Lucianna Calia for help in the experiments reported in Figure 1—figure supplement 11 and in Figure 5—figure supplement 2, respectively. We are in debt with Nacho Molina and Davide Mazza for suggestions and critical reading of the manuscript. We also thank the Millipore service for assistance in microfluidics set-up. We thank the reviewing team for their help and constructive comments.

Reviewing Editor

  1. Suzanne Gaudet, Reviewing Editor, Dana-Farber Cancer Institute, Harvard University, United States

Publication history

  1. Received: June 1, 2015
  2. Accepted: January 13, 2016
  3. Accepted Manuscript published: January 14, 2016 (version 1)
  4. Version of Record published: March 10, 2016 (version 2)

Copyright

© 2016, Zambrano et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,377
    Page views
  • 701
    Downloads
  • 15
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Cell Biology
    2. Genes and Chromosomes
    Wahid A Mulla et al.
    Research Article Updated