Understanding the dynamics and variability of protein circuitry requires accurate measurements in living cells as well as theoretical models. To address this, we employed one of the best‐studied protein circuits in human cells, the negative feedback loop between the tumor suppressor p53 and the oncogene Mdm2. We measured the dynamics of fluorescently tagged p53 and Mdm2 over several days in individual living cells. We found that isogenic cells in the same environment behaved in highly variable ways following DNA‐damaging gamma irradiation: some cells showed undamped oscillations for at least 3 days (more than 10 peaks). The amplitude of the oscillations was much more variable than the period. Sister cells continued to oscillate in a correlated way after cell division, but lost correlation after about 11 h on average. Other cells showed low‐frequency fluctuations that did not resemble oscillations. We also analyzed different families of mathematical models of the system, including a novel checkpoint mechanism. The models point to the possible source of the variability in the oscillations: low‐frequency noise in protein production rates, rather than noise in other parameters such as degradation rates. This study provides a view of the extensive variability of the behavior of a protein circuit in living human cells, both from cell to cell and in the same cell over time.
A major challenge of systems biology is to understand the dynamics of regulatory circuits (Hartwell et al, 1999; Becskei and Serrano, 2000). In order to understand such circuits, it is important to study them at the level of individual cells, rather than in averages of cell populations (Ferrell and Machleder, 1998). In this study, we investigated the dynamics of one of the most studied protein regulatory circuits in human cells—the p53–Mdm2 negative feedback loop (Lev Bar‐Or et al, 2000; Vogelstein et al, 2000; Oren et al, 2002; Oren, 2003; Harris and Levine, 2005) in individual living cells. p53 is a tumor suppressor that upregulates a variety of genes that are involved in DNA repair and cell proliferation (Vogelstein et al, 2000; Oren, 2003). Interestingly, one of p53's downstream targets is Mdm2, which inhibits p53 activity and targets it for degradation. This negative feedback loop was previously shown, in a study that followed cells over 16 h, to exhibit one to two discrete undamped peaks following DNA‐damaging radiation (Lahav et al, 2004). Such detailed information is lost in population studies, as the peaks in different cells are averaged together, giving an appearance of damped oscillations.
Here we extended our previous study on individual living cells for much longer periods of up to 3 days. In order to trigger the p53–Mdm2 response, we used gamma irradiation that causes double‐stranded DNA breaks. We found that the p53/Mdm2 oscillations are continuous for at least 60 h: many cells show more than 10 oscillation peaks (Figure 2). We found that the oscillation pattern was highly variable between cells that were genetically identical. In addition to cells that oscillated, other cells showed a dynamic fluctuation of protein levels that did not resemble sustained oscillations. The prolonged experiments indicate that the fraction of cells that show continuous oscillations increases with the irradiation dose: in our previous, shorter experiments, this effect could not be seen owing to the small number of peaks.
During the extended duration time‐lapse microscopy of live cells, we observed cell divisions and could address the dynamics of the two daughter cells. We found that sister cells continued to oscillate in a correlated way after cell division, but lost correlation after about 11 h on average. We also found that unirradiated cells show slow fluctuations in p53 and Mdm2 levels, with frequencies on the order of 11 h.
Finally, we studied several families of mathematical models in order to understand the variability in the dynamics observed in our experiments (Figure 6). As many of the details of the system are unknown, we studied six different model designs based on known interactions and suggested models of the p53 system (Tiana et al, 2002; Tyson et al, 2003; Ciliberto et al, 2005; Ma et al, 2005). The main conclusions of the modeling study are as follows:
Different models can describe deterministic oscillations, and do so with very similar parameters: ‘consensus parameters’. These parameters, such as p53 and Mdm2 lifetimes, can be directly tested experimentally.
We went beyond deterministic models, and performed an analysis of the p53 system with stochasticity in the internal parameters. We find that the best way to get the observed noise in amplitude but relative precision in the period is to add noise to the production rates, rather than degradation rates. This might point to the internal source of the noise in this system.
We performed an analysis of internal noise that has slow variation (rather than rapidly varying ‘white noise’ usually employed in most models of biological noise). We find that such slow noise is essential to explain the behavior observed: noise that is too fast or too slow is averaged out and cannot give significant variations in the oscillations. It seems that only noise with frequencies near the resonance frequency of the oscillator can result in the observed amplitude variations.
We find that a new checkpoint model captures all of the features of the variability we observe, including low correlation between p53 and Mdm2 peak amplitudes. This model predicts that a factor downstream of p53 inhibits a factor upstream of p53, such as phosphorylated ATM.
In summary, this is one of the first studies of a protein circuit as it functions in individual living human cells. We find unexpected variability, prolonged oscillations after input signals, and unsynchronized fluctuations in the absence of input signal. All of these could not be seen in experiments that average over cell populations. Models help to pinpoint the intracellular source of these variations.
Irradiation results in prolonged oscillation in the p53 system, lasting for at least 3 days with a relatively defined period of about 6 hours but noisy amplitude.
Some cells show unsynchronized, low‐frequency (∼11 hour) fluctuations in p53 and Mdm2 levels in irradiated cells as well as in non‐irradiated cells.
In order for models to capture the observed variability, they must include slowly‐varying fluctuations in protein production rates, rather than noise in other parameters.
A goal of systems biology is to understand the design principles that govern dynamics of protein regulatory circuits (Hartwell et al, 1999). It is especially important to investigate network motifs, regulatory patterns that recur in various biological networks (Milo et al, 2002; Alon, 2003, 2006). Understanding the dynamical features of a specific network motif may help us to understand diverse biological systems in which this motif appears (Lee et al, 2002; Shen‐Orr et al, 2002; Mangan and Alon, 2003; Kalir and Alon, 2004; Odom et al, 2004; Boyer et al, 2005; Ma'ayan et al, 2005; Mangan et al, 2006). For this purpose, it is important to study the best‐characterized systems using dynamic measurements in living cells.
To understand protein circuits, it is important to study the impact of the stochastic nature of biological reactions on the behavior of the circuit (Novick and Weiner, 1957; Spudich and Koshland, 1976; McAdams and Arkin, 1997, 1999; Becskei and Serrano, 2000; Thattai and van Oudenaarden, 2001; Elowitz et al, 2002; Hasty et al, 2002; Ozbudak et al, 2002; Blake et al, 2003; Isaacs et al, 2003; Paulsson, 2004; Raser and O'Shea, 2004, 2005; Becskei et al, 2005; Blake and Collins, 2005; Colman‐Lerner et al, 2005; Golding et al, 2005; Kaern et al, 2005; Sachs et al, 2005; Weinberger et al, 2005; Volfson et al, 2006). For this purpose, it is essential to study individual cells and to measure the cell–cell variations in the biological response, rather than averaging over cell populations. Most studies of stochastic behavior to date have been in microorganisms. It would therefore be of interest to measure the behavior and variability of a protein circuit over long times in individual human cells.
Here, we study the dynamics and variability of one of the network motifs that recurs across organisms: a negative feedback loop, which is composed of interactions on two different timescales—a slow positive transcriptional arm and a fast negative protein–protein interaction arm (Lahav et al, 2004; Yeger‐Lotem et al, 2004; Ma'ayan et al, 2005). We study this network motif within one of the best‐characterized systems in human cells, the negative feedback loop between p53 and Mdm2 (Kubbutat and Vousden, 1998; Prives, 1998; Larkin and Jackson, 1999; Prives and Hall, 1999; Vogelstein et al, 2000; Ryan et al, 2001; Vousden and Lu, 2002; Oren, 2003; Meek, 2004; Bond et al, 2005; Harris and Levine, 2005).
In the p53 system, p53 transcriptionally activates mdm2. Mdm2, in turn, negatively regulates p53 by both inhibiting its activity as a transcription factor and by enhancing its degradation rate (Barak et al, 1993; Wu et al, 1993; Haupt et al, 1997; Kubbutat et al, 1997; Piette et al, 1997; Momand et al, 2000). The concentration of p53 increases in response to stress signals, such as DNA damage. The main mechanism for this increase is stabilization of p53 due to reduced interaction with Mdm2. Following stress signals, p53 activates transcription of several hundred genes that are involved in growth arrest, apoptosis, senescence, and DNA repair. It is important to note that many additional proteins interact with p53 and Mdm2, so that the negative feedback loop motif is embedded inside a network of additional interactions, many of which are not fully characterized (Harris and Levine, 2005).
Models of negative feedback loops, such as between p53 and Mdm2, suggest that they can generate an oscillatory behavior with a time delay between the two proteins (Lev Bar‐Or et al, 2000; Mihalas et al, 2000; Hoffmann et al, 2002; Tiana et al, 2002; Michael and Oren, 2003; Monk, 2003; Tyson et al, 2003; Nelson et al, 2004; Tyson, 2004; Ciliberto et al, 2005; Ma et al, 2005). For different parameters of the feedback loop, the dynamics can show either a monotonic response, damped oscillations, or undamped (sustained) oscillations in which each peak has the same amplitude as the previous peak (Lahav, 2004). The stronger the interactions between the proteins, the more oscillatory the dynamics. Other parameters, such as high basal degradation rates of the proteins, tend to damp out the oscillations. Most models of the p53 network used deterministic equations, and thus did not consider the cell–cell variability in the dynamics.
Experimental studies have shown that p53 and Mdm2 undergo oscillatory behavior following DNA damage caused by gamma irradiation. These oscillations appeared to be damped in assays that measure averages over population of cells (Lev Bar‐Or et al, 2000). In a previous study, we developed a system for following p53 and Mdm2 dynamics in individual living cells. This system used an MCF7 breast cancer cell line stably transfected with p53 fused to cyan fluorescent protein (CFP), and Mdm2 fused to yellow fluorescent protein (YFP). The p53‐CFP fusion protein was active in causing apoptosis and transactivating downstream targets. The concentrations and dynamics of both fluorescently tagged proteins were found in Western blots to reliably reproduce the concentration and dynamics of the endogenous p53 and Mdm2 expressed by these cells (Lahav et al, 2004). In our previous study, individual cell measurements using fluorescent microscopy were limited to 16 h following gamma irradiation. During this 16 h period, we observed up to two undamped peaks of p53‐CFP concentration following gamma irradiation (Lahav et al, 2004). The peak amplitude and timing did not depend on the dose of irradiation. The mean number of peaks appeared to increase with irradiation dose, in the sense that the probability for two peaks in the 16 h experiment increased with dose, whereas the probability for no oscillation peaks decreased with dose.
Here, we experimentally investigate p53 and Mdm2 dynamics in individual living breast cancer (MCF7) cells for much longer times than in our previous study. In a large fraction of cells, we find sustained undamped oscillations of p53‐CFP and Mdm2‐YFP, which lasted for at least ∼3 days following gamma irradiation. We also extend our previous study by examining the noise in the response. We find that the oscillation pattern was highly variable between isogenic cells, but this variation had distinct properties: the oscillation amplitudes fluctuated widely, yet the oscillation frequency was much less variable. In addition to cells that oscillated, other cells showed a dynamic fluctuation of protein levels that did not resemble sustained oscillations. The prolonged experiments indicate that the fraction of oscillating cells increases with irradiation dose.
We also present a theoretical analysis of the negative feedback loop in the p53 system. We extend previous theoretical studies by investigating several families of models, and by studying the effect of stochasticity in the model reactions. We find that several distinct model families can capture the experimentally observed oscillations, and suggest ‘consensus parameters’ for in vivo degradation and production rates in this system. To capture the variability in the dynamics, we find that one must explicitly add long‐wavelength noise to different model parameters. The analysis indicates that the observed characteristic variation in the oscillations stems from fluctuations in the protein production rates, rather than from fluctuations in other parameters. Essentially, the negative feedback loop amplifies slowly varying noise in the protein production rates at frequencies near the resonance frequency of the feedback loop.
Prolonged time‐lapse movies show undamped oscillations over days
We used an MCF7 cell line clone stably transfected with p53‐CFP and Mdm2‐YFP (Lahav et al, 2004). The results we describe were obtained from isogenic cells, grown from a single parental cell (Lahav et al, 2004). Western blots indicate that the concentration of exogenous p53‐CFP and Mdm2‐YFP protein in our cell line is comparable to the endogenous p53 and Mdm2 protein concentration (Lahav et al, 2004). Hence, these proteins are not strongly overexpressed in the present system.
We obtained time‐lapse fluorescence microscopy movies of these cells for extended periods of time after exposure to gamma irradiation (Figure 1A and movie SM1 in Supplementary information). Overall, we collected time courses from over 1000 individual cells in different experiments with different doses of gamma irradiation. Most of the time‐lapse movies were performed in an incubator environment with controlled humidity, temperature, and CO2, providing conditions that allow growth over several days. Every 10–20 min, images of the cells in fluorescence and phase illumination were captured. Cells divided vigorously in the movies without gamma irradiation for at least 3 days. Gamma irradiation caused cells to enter growth arrest. We found that p53‐CFP and Mdm2‐YFP fluorescence was brightly visible when the proteins were in the cell nucleus (Figure 1A and Supplementary Figure S1). For each cell, we obtained a time‐dependent signal equal to the mean fluorescence intensity of p53‐CFP and Mdm2‐YFP in the nucleus (Figures 1 and 2 and Supplementary Figure S2).
Nuclear levels of p53‐CFP and Mdm2‐YFP were found to oscillate continuously following gamma irradiation in a large number of cells. In most of these cells, these oscillations lasted for the entire movie duration (Figures 1 and 2).
We analyzed the characteristic oscillation frequency in each cell using Fourier analysis (Supplementary Figure S3) and pitch detection, a method commonly used for determining principal frequencies in noisy non‐stationary signals in the context of speech recognition (see Materials and methods). In long movies, we found that about 60% of the cells exposed to a pulse of 10 Gy gamma irradiation showed sustained Mdm2‐YFP oscillations, with a period of 5.5±1.5 h (Figure 3).
It is important to note that a significant fraction of the cells (about 40% in 10 Gy) showed Mdm2‐YFP dynamics that did not resemble sustained oscillations (Figure 2 and Supplementary Figure S2). These cells showed either no response or slowly varying fluctuations (Figure 2A and B, bottom panels). In a few cells, the oscillations stopped or changed frequency after 1–2 days.
The onset of oscillations in different cells was synchronized to the DNA damage signal. Cells gradually lost synchrony with each other owing to the variations in oscillation frequencies (Figure 2D and Supplementary Figure S2E). In oscillating cells, Mdm2‐YFP peaks followed p53‐CFP peaks at a delay of 2±0.5 h on average (Figures 1 and 4C).
We evaluated the amplitude and width of each peak in each oscillating cell, and calculated the average of these properties. The average amplitude of the oscillations did not appear to change significantly over time (Figure 4A). Similarly, the mean peak width did not change considerably throughout the movies (Figure 4B). In this sense, the oscillations can be described as undamped.
Peak amplitude is highly variable, whereas peak timing is more precise
The dynamics of cells from a clone in the same field of view showed significant cell–cell differences. These differences were seen between different cells, and also between different peaks in the same cell. We examined the variability between peaks in the oscillations (e.g. Figure 1B). We found that the amplitudes of the individual peaks varied with a coefficient of variance (standard deviation divided by mean) of about 70% (Figure 4D). The amplitudes of Mdm2‐YFP peaks were not correlated to the amplitude of the preceding or the subsequent p53‐CFP peaks (correlation coefficient of 0±0.2). In some cases, Mdm2‐YFP peaks occurred without a detectable preceding p53‐CFP peak (Supplementary Figure S4).
In contrast to the large variability in amplitude, the peak width and p53–Mdm2 delay of individual peaks were more constant and varied by only about 30% (Figure 4E and F). The variation in the oscillation period for each cell (change of pitch value along the oscillation signal) was less than 20% in most oscillating cells.
Correlation between sister cells is lost within half a generation
To further study the variability in the dynamics of each cell, we examined cells that underwent cell division during the movie. In the first 40 h of the movies, these comprised ∼75% of the non‐irradiated cells, ∼65% of the cells following 0.3 Gy, ∼50% of the cells following 5 Gy, and ∼10% of the cells following 10 Gy of gamma irradiation. Six typical examples are shown in 3, 4Figure 5A.
We analyzed over 100 sister‐cell pairs following cell division (see Materials and methods). We found that after division, the dynamics of Mdm2‐YFP were correlated between sister cells for a few hours. This correlation was reduced by 50% within about 11±5 h on average (Figure 5B). No significant correlation of the dynamics with the cell cycle was observed in these cells.
Some cells show non‐oscillatory fluctuations
We found that a fraction of both irradiated cells and non‐irradiated cells showed Mdm2‐YFP signals that had slowly varying fluctuations that did not resemble oscillations (Figure 2 and Supplementary Figure S2). The fluctuations were rather slow, with a typical timescales of 8–12 h, as determined by Fourier analysis and pitch‐detection methods (Figure 3). The fluctuations showed at most 2–3 such peaks rather than sustained oscillations (Figure 2A–C and Supplementary Figure S2). Similar fluctuations were also observed in p53‐CFP (data not shown). Control cells expressing a YFP fusion protein (YFP fused to histone H2AZ), showed no such fluctuations (data not shown).
Fraction of cells with ∼5.5‐h oscillations increases with gamma dose
We also measured the dynamics under different doses of irradiation. When gamma irradiation was applied at 0.3 or 5 Gy, a fraction of the cells displayed oscillations with a period of about 5.5 h, similar to those in the 10 Gy experiment. We used pitch detection to estimate the fraction of cells whose characteristic period is 4–7 h. We found that the fraction of cells that perform Mdm2‐YFP oscillations increases with gamma dosage (Figure 3). For all irradiation doses, the oscillations in these cells typically showed many peaks and were undamped. The mean amplitude and period of oscillations in individual cells did not appear to significantly depend on irradiation level 5(Supplementary Figure S5).
Several families of models can generate the observed oscillations and show consensus biochemical parameter values
We considered several mathematical models of the p53–Mdm2 feedback loop. Since current knowledge of the system is incomplete, we analyzed the simplest possible models, aiming to understand the general properties of model families. Our motivation was to find a simple model that could capture the characteristics of the undamped and noisy oscillations that were found in many of the cells.
We examined six model families (Figure 6A and Table I). All of the models include the negative feedback loop in which p53, denoted by x, transcriptionally activates Mdm2 denoted by y. Active Mdm2 increases the degradation rate of p53.
Three of the models are delay oscillators (I, III and IV) (Mihalas et al, 2000; Goldbeter, 2002; Tiana et al, 2002; Monk, 2003). The models differ in mathematical details that describe the delay between x and y and the effects of y on x. Model I includes an Mdm2 precursor, denoted by y0, representing, for example, Mdm2 mRNA, and the action of y on x is described by first‐order kinetics in both x and y. In model IV, the action of y on x is nonlinear, and described by a saturating Michaelis–Menten function. In model III, the Mdm2 precursor y0 is replaced by a stiff delay term, which makes the production rate of Mdm2 depend directly on the concentration of p53 at an earlier time. A recent model by Ma et al (2005) and Wagner et al (2005) combines features of models III and IV.
In addition to the three delay oscillators, we also considered two relaxation oscillators (II and V) (Wilhelm and Heinrich, 1995; Murray, 2003; Pomerening et al, 2003; Tyson et al, 2003; Ciliberto et al, 2005). In these models, the negative feedback loop is supplemented by a positive feedback loop on p53. This positive feedback loop might represent in a simplified manner the action of additional p53 system components, which have a total upregulating effect on p53 (Harris and Levine, 2005). This type of model was recently studied by Ciliberto et al (2005). We considered both linear positive regulation (model V) and nonlinear regulation based on a saturating function (model II).
These models (I–V), although differing in detail, rely on a single negative feedback loop. The last model (VI) is a novel checkpoint mechanism, which uses two negative feedback loops, one direct feedback and one longer loop that impinges on an upstream regulator of p53. In this model, a protein downstream of p53 inhibits a signaling protein that is upstream of p53 (see more details in Supplementary information; Banin et al, 1998). For simplicity, this inhibitor is modeled by y, but it need not be Mdm2 and can also represent a different protein with similar dynamics. This model predicts that upstream elements (e.g. phosphorylated ATM) also undergo oscillatory dynamics (see appendix in Supplementary information and Supplementary Figure S9). The model was inspired by the observation that an upstream regulator of p53, namely phosphorylated ATM (Bakkenist and Kastan, 2003) that responds to double‐stranded DNA breaks (DSBs), shows a pulse of activity after application of a radiomimetic drug (NCS) in a set of Western blot experiments measuring protein dynamics for 6 h after damage (Banin et al, 1998; Stommel and Wahl, 2004).
We numerically solved all six models for a wide range of parameters. We selected parameter values (Table II) that best reproduce an effective averaged individual cell measurements of nuclear p53‐CFP and Mdm2‐YFP following gamma irradiation, that is, oscillations that do not dampen out considerably and have constant inter‐peak timing (Supplementary Figure S6). Note that this effective individual cell dynamics is different from the population‐average dynamics, which is a damped oscillation (Supplementary Figure S2E).
Model I cannot produce multiple oscillations similar to those experimentally observed. The oscillations in models II and III are very sensitive to parameters. Small changes in some of the parameters listed in Table II cause these models to show strongly damped oscillations (Supplementary Figure S8). Such sensitive (non‐robust) circuits might not be expected to function properly in the noisy cellular context (Savageau, 1976; Barkai and Leibler, 1997; Alon et al, 1999; Eldar et al, 2002; Kitano, 2004).
In contrast, models IV–VI could generate sustained or weakly damped oscillations (Figure 6B and Supplementary Figure S6) over a broad range of parameters (Supplementary Figure S8). Importantly, most of the parameters shared by all three models showed very similar best‐fit values. Thus, these models may provide estimates of the effective biochemical parameters such as production rates and degradation times of p53 and Mdm2. The ‘consensus’ values of the parameters are shown in Table II. In all three models, the Mdm2 degradation rate was about αy≈1 h−1, the time for Mdm2 maturation was about 1/α0≈τ≈1 h, and the Mdm2‐independant degradation rate of p53, αx, was negligible. The models also agreed on the values βx and βy, the rates of p53 and Mdm2 production.
Increasing the Mdm2 lifetime or its maturation time led, according to the models, to a lower natural frequency and to pulses with longer time periods (Supplementary Figure S8). This might help explain the low‐frequency fluctuations observed with no gamma irradiation (Figure 2C), because Mdm2 lifetime is longer in the absence of DNA damage than in its presence (Stommel and Wahl, 2004).
The observed noise in oscillation amplitude is captured by low‐frequency fluctuations in the protein production rates
Deterministic simulations cannot capture the variability in the oscillation amplitudes observed in the cells (Figures 1, 2, and 4). We therefore added internal stochasticity to the equations. We found that the characteristic variability observed in our experiments, where amplitude varies more strongly than frequency, could best be captured by varying the protein production rates. Production rate variations change amplitudes, but do not significantly affect the oscillation period (Supplementary Figure S8). In contrast, we find that variations in other parameters, such as degradation rates, generally lead to variations in both amplitude and period.
To describe stochasticity in protein production rates, we used multiplicative noise in the protein production terms. Most previous theoretical analyses of noise employed white noise, which is rapidly fluctuating (McAdams and Arkin, 1997; Thattai and van Oudenaarden, 2001; Paulsson, 2004; Kaern et al, 2005; Ramanathan and Swain, 2005). We used, in addition to white noise, noise with different characteristic correlation times, including noise that varies on slow timescales. This was inspired by the recent observation that protein production rates vary significantly between individual bacterial cells, and that this variation has long autocorrelation times on the order of a cell generation (Rosenfeld et al, 2005).
We first used high‐frequency noise (similar to white noise), in which the production rates varied with a correlation time on the order of minutes to an hour (Figure 6C). This may represent intrinsic noise due to stochastic transcription and translation (McAdams and Arkin, 1997; Thattai and van Oudenaarden, 2001; Elowitz et al, 2002; Ozbudak et al, 2002; Blake et al, 2003; Isaacs et al, 2003; Paulsson, 2004; Becskei et al, 2005; Colman‐Lerner et al, 2005; Golding et al, 2005; Kaern et al, 2005; Volfson et al, 2006). The stronger the noise in production rates, the higher the resulting fluctuations in the dynamics. However, even very strong high‐frequency noise (STD of 50% in production rates) resulted in only mild variations in the oscillation amplitudes in all six models. These variations were significantly lower than the variations observed in individual cell measurements in the present experiments.
We next introduced low‐frequency noise, with a timescale of several hours (e.g. 12.5 h; Figure 6D). We found that under such noise, the amplitudes vary far more strongly than under high‐frequency noise. The extent of the variability was similar to that experimentally observed, with strong amplitude variations and smaller variations in the period of the oscillations.2
Finally, we found that very low‐frequency noise (e.g. 50 h; Figure 6E) does not produce strong variability in the oscillations. It appears that the oscillators can only amplify the frequency component of the noise close to their natural resonant frequency (about 6 h). We find that the variability is maximal at noise frequencies of about twice the natural frequency of the oscillator (Figure 6F), such that consecutive peaks are oppositely affected by the noise.
The models with low‐frequency noise in the production rates showed qualitatively similar dynamics to those found in the experiments, including occasional loss of a peak. Only model VI was able to reproduce our observations that p53 and Mdm2 peak amplitudes had only a weak correlation. Other models had a strong coupling in the variations of the peaks of these two proteins.
The present study examined p53 and Mdm2 dynamics in individual cells from a clone. We found undamped oscillations with more than 10 consecutive peaks, lasting for at least three days following DNA damage. The dynamics showed striking cell–cell variability. A fraction of the cells showed either no response or a slowly fluctuating signal. The cells that performed oscillations displayed large variation in peak amplitude, and smaller variations in the oscillation period. Models point to the source of the noise in the oscillations: low‐frequency fluctuations in protein production rates.
The oscillations following DNA‐damaging gamma irradiation had a period of about 5.5 h, and were synchronized to the gamma irradiation pulse. The number of oscillating cells increased with gamma dose, reaching about 60% of the cells following 10 Gy. Some cells divided during the movies. Divisions allowed one to follow the passage of information across the cell division event. We found that the oscillations continued in the same phase after division, suggesting that the information in the system is transferred to the daughter cells. However, correlation between daughter cells was lost after about 11 h. This loss of correlation indicates the timescale on which prediction of the cell state in this system can be made based on the cell state in the past.
The oscillations after DNA damage had a distinct noise characteristic. Their amplitudes varied from peak to peak by about 70%. In contrast to the large amplitude variation, the oscillation period was less noisy, and had a variability of only 20%. Similar features are seen in other biological oscillators. For example, the cell‐autonomous circadian clock in cyanobacteria and in fibroblasts shows larger amplitude variations than timing variations in experiments and in models (Barkai and Leibler, 2000; Vilar et al, 2002; Mihalcescu et al, 2004; Nagoshi et al, 2004). Precise period and variable amplitude may characterize other biological oscillators.
Although the timing is relatively precise, and the oscillations are initially synchronized to the gamma irradiation signal, the variation in timing causes peaks to eventually go out of phase. Therefore, the p53 and Mdm2 dynamics appear as damped oscillations in assays that average over cell populations, such as immunoblots. This averaging effect was also seen by averaging over the present individual cell dynamics, showing damped oscillations with 2–3 discernable peaks (Supplementary Figure S2E).
It is interesting to compare the present results with our previous study that followed cells over only 16 h (Lahav et al, 2004). In that study, cells showed either zero, one, or two peaks of p53 in the 16 h period. The fraction of cells with two peaks increased with gamma irradiation. It seemed therefore that the number of peaks depended on the gamma dose. The present study, which followed cells over a much longer time, suggests that oscillations in most cells are in fact long lasting, and that most oscillating cells show numerous peaks following damage. We found that the fraction of oscillating cells (with a 4–7 h period in Mdm2‐YFP levels) increases with gamma dose. The previous 16 h movies registered some cells with one pulse, whereas the present study indicates that such cells can often show additional pulses after a delay (Supplementary Figure S4). This emphasizes the importance of extended measurements for dynamical systems with slow timescales.
How are the oscillations produced? Instead of analyzing a single model, the limited state of current knowledge of the system makes it appropriate to study several families of models, to ask about the general properties of the dynamics. We performed a theoretical analysis of several model families. Most models were able to produce oscillations. The models suggest that the noise in the oscillations is owing to stochasticity in the protein production rates, rather than in other parameters such as degradation rates. Furthermore, the observed oscillations suggest that the noise in protein production rate has a slowly varying component, with a correlation time of 10–20 h. Internal noise that is too fast or too slow cannot explain the observed variability. The negative feedback loop, which is a natural oscillator, amplifies the frequency component of the noise in the vicinity of its natural frequency, resulting in the observed variability.
The present results were obtained in a clonal population of a human, MCF7 cell line, stably expressing fluorescent fusions of p53 and Mdm2. Endogenous p53 and Mdm2 oscillations were found in cell averages also in MCF7 cells that do not express ectopic fusion proteins (Lev Bar‐Or et al, 2000). These cancer cells might be deficient in some aspects of p53 regulation (Vojtesek and Lane, 1993) and downstream apoptotic responses (Janicke et al, 1998). It would therefore be important to study other cell types. For example, Western blots performed over several hours after DNA damage showed a peak of p53, Mdm2, and p21 expression in several cell lines including WS1 human primary skin fibroblasts (Stommel and Wahl, 2004), human glioblastoma cells (Ohnishi et al, 1999), and HCT116 human colon cancer cells (Chen et al, 2005). This raises the possibility that oscillations may occur also in these cell lines. It would be important to extend the present individual cell experiments to other cell types, and to try to monitor DNA damage in parallel to the dynamics of the p53 system.
Perhaps the most intriguing question raised by these observations is the biological function of the undamped oscillations, assuming that they also occur in normal cells with endogenous p53 and Mdm2. One clue is that undamped oscillations are also found in other stress‐response systems. Tightly regulated oscillations with variable amplitude and precise timing were recently observed in the SOS DNA‐damage response in Escherichia coli (Friedman et al, 2005). Highly variable nuclear‐cytoplasmic oscillations were found in NF‐κB system (Hoffmann et al, 2002; Nelson et al, 2004). Both NF‐κB and the SOS regulator LexA are involved in a negative feedback loop motif similar to that of p53–Mdm2. As in the p53 system, these loops are embedded in many additional interactions. The presence of oscillations in the systems mentioned above may suggest that oscillations play a general role in stress or damage response.
The present study demonstrated prolonged undamped oscillations in the p53–Mdm2 system following gamma irradiation. Significant cell–cell variability was observed in the amplitude but not period of the oscillations. Some of the cells had slow fluctuations that do not resemble oscillations; the fraction of oscillating cells increased with irradiation dose but the oscillation amplitude did not. Modeling suggests that the noise in the oscillations reflects slow internal noise in protein production rates. The present approach that combines long‐term dynamic experiments in individual cells and theoretical analysis of families of models may help to understand oscillations and cell–cell variability in other regulatory systems.
Materials and methods
Cell line and constructs
We used MCF7, human breast cancer epithelial cells, U280, stably transfected with pU265 and pU293 as described (Lahav et al, 2004). In pU265, ECFP from pECFP‐C1 (Clontech) was subcloned after the last codon of p53 cDNA, under the mouse Metallothionein‐1 promoter (MTΔ156) (Brinster et al, 1982). This promoter provides a basal and constant level of transcription of p53‐CFP. A basal promoter for p53‐CFP was chosen because p53 is thought to be primarily regulated at the protein level and not at the transcriptional level (Michael and Oren, 2003). Control experiments with CFP expressed from this promoter showed constant expression with no oscillations. In pU293, the hMDM2 promoter was cloned by PCR using genomic DNA as a template, creating a 3.5 kb fragment upstream of the ATG site in exon 3, including P1 and P2 (Oliner et al, 1992). This promoter was subcloned into pEYFP‐1 (Clontech) (Lahav, 2004).
Cells were maintained at 37°C in 96‐well plates or in 2 mm optical plates (Nunc) in RPMI 1640 medium containing 10% fetal calf serum (Sigma). At 1–2 h before observation in the microscope, medium was changed to RPMI 1640 medium containing 3% fetal calf serum, HEPES, and 2 mM l‐glutamine, lacking riboflavin and phenol red (Beit Haemek, Biological Industries), in order to reduce background fluorescence. Cells were then exposed to the appropriate dose of gamma irradiation (60Co, 1.8 Gy min−1). The number of DSBs has been found to be linear in gamma dose, with an average of about 30 DSB per Gy per cell (Bonner, 2003). Cells were viewed with three types of inverted fluorescence microscope systems denoted by MS.I, MS.II, and MS.III. MS.I: Olympus IX70 with a Photometrics Quantix 57 cooled back‐illuminated CCD camera, in a 37°C incubator, using bright‐field, CFP and YFP exposures, every 20 min, with a mercury lamp. MS.II: Leica DMIRE2 with a Hamamatsu ORCA‐ER cooled back‐illuminated CCD camera, in a 37°C incubator with humidity and CO2 control, using phase‐contrast and YFP exposures only, every 10 min, with a mercury lamp. MS.III: Nikon TE2000E2 with a Hamamatsu ORCA‐ER cooled back‐illuminated CCD camera, in a 37°C incubator with humidity and CO2 control, using phase‐contrast, YFP and CFP exposures, every 20 min, with a xenon lamp.
CFP filter set: excitation 436/20 nm, dichroic beam splitter 455 nm, emission 480/40 nm. YFP filter set: excitation 500/20 nm, a dichroic beam splitter 515 nm, emission 535/30 nm.
The mean cell generation time was about 20 h in the CO2 incubated microscope without gamma irradiation. We find that movies using CFP and YFP illumination over 3 days did not visibly affect the cell morphology or generation time.
Cell tracking and fluorescence quantification
Cell images captured in MS.I (Figure 1 and Supplementary Figure S4): Relative fluorescence analysis and background subtraction was carried out using custom written Matlab software (Mathworks Inc.). The location of each cell nucleus was marked manually in each frame, using a custom written graphical user interface in Matlab. Independent tracking by four different researchers showed that this manual step contributed <5% errors. Background fluorescence was measured at manually marked locations with no cells, and subtracted from the nuclear fluorescence. Mean fluorescence intensity of pixels in the nucleus was measured. Cellular autofluorescence of wild‐type MCF7 cells without the CFP or YFP genes gave consistent and low values with a mean of 25 CFP fluorescence units per pixel and 1 YFP fluorescence unit per pixel, with a coefficient of variation of ∼30%. In these units, average peak amplitude (range from minimum to maximum) was ∼45 CFP fluorescence units (for p53‐CFP) and ∼8 YFP fluorescence units (for Mdm2‐YFP).
Cell images captured with MS.II (Figures 2 and 5, and Supplementary Figure S2) and MS.III: Relative fluorescence analysis and background subtraction was carried out using custom written Matlab software (Mathworks Inc.). Nuclei identification and tracking was performed using MetaMorph™ software, and was manually controlled for accuracy. Comparison of sister cells (that were separately tracked) on frames before cell division shows that the identification and tracking contributes ∼2% errors. Background was automatically subtracted. Mean fluorescence intensity of pixels in the nucleus was measured. Autofluorescence (in the YFP channel) was negligible. The inhomogeneity of the illumination was measured using a solution of purified GFP (BD Biosciences Clontech, Palo Alto, CA) before and after every movie and automatically corrected using custom written Matlab software (Mathworks Inc.). Bleaching effects were corrected using an empirical fit of the mean nuclear fluorescence levels to a decaying exponent with an offset. Independent controls, in which H1299 cells with constitutive nuclear YFP expression were imaged, indicate that measurement errors and fluctuations in this system are on the order of a few percent.
Statistical analysis of pulse properties
The p53‐CFP and Mdm2‐YFP data from MS.I were analyzed in the time domain. In the dynamic curve of each cell, separate pulses of expression were manually marked using custom written software (Matlab). The separate pulses were identified using criteria based on pulse magnitude and signal‐to‐noise ratio. The baseline was subtracted for each pulse separately, to correct for slowly varying noise in the fluorescence quantification, which may originate from slowly varying autofluorescence of cells. For the comparison of the pulse properties, the time domain was divided into segments of length 300–400 min, and each pulse was independently assigned an ordinal number according to the time segment when it occurs. Average properties (and standard errors) were then calculated for all the pulses that occur at a given time interval (Figure 4 and Supplementary Figure S5).
The intensity signals from each cell, as obtained from the microscope, were analyzed to detect the oscillation period (1/frequency). We used a standard method for the detection of pitch, used in speech and music processing (Rabiner and Schafer, 1978). Pitch can be considered as the basic frequency of oscillations. Each signal was divided into segments of 128 samples, with a sliding window, which was moved at increments of eight samples. For each window, the autocorrelation of the windowed segment was computed, and normalized so that the autocorrelations at zero lag are identically 1. The first peak of the autocorrelation function was detected and identified as the pitch period of this window if its autocorrelation value was higher than 0.2. The sliding window method enables tracing temporal changes in the oscillation period. To detect the most prominent pitch period for each cell, we binned separate segment periods into 10 bins and selected the most common period.
At each time point, we ranked all the cells in a movie from lowest to highest nuclear Mdm2‐YFP fluorescence level, and normalized the ranking to the range 0–1. For a random pair of cells, the absolute difference in ranking is equal to Dr=1/3 on average. For each pair of sister cells (after division), we measured the absolute difference in rank between the two sister cells. We calculated the average for all sister pairs as a function of time after division, for a total of 112 sister‐cell pairs after 0.3, 5, and 10 Gy of gamma irradiation. This average difference, D(t) was found to increase over time from D(t=0)∼0.05 (equivalent to the minimal rank difference in a movie with 20 cells) to D(t>30 h)∼0.3, close to the population average between unrelated cells. In Figure 5, we plot the normalized sister‐pair difference, D′(t)=(Dr−D(t))/(Dr−D(0)). Similar results (half correlation time of 6–16 h) were found with different measures of average sister‐cell rank differences, such as root‐mean‐square difference, and with different subsets of cell (such as only those exposed to 0.3 or 5 Gy).
Numerical integration and optimization were carried out using Matlab software.
We thank Michael B Elowitz for assistance and contributions, Arnold J Levine and Moshe Oren for encouragement discussions and suggestions, Gareth L Bond, Gustavo A Stolovitzky, Lan Ma, John M Wagner and all members of our labs for comments and discussions. We acknowledge support grants from the European Commission (COMBIO: LSHGCT‐2004‐503568) and from the Kahn Fund for Systems Biology at the Weizmann Institute of Science.
Supplementary Figures and Information
- Copyright © 2006 EMBO and Nature Publishing Group