Noise in the Vertebrate Segmentation Clock Is Boosted by Time Delays but Tamed by Notch Signaling

SUMMARY Taming cell-to-cell variability in gene expression is critical for precise pattern formation during embryonic development. To investigate the source and buffering mechanism of expression variability, we studied a biological clock, the vertebrate segmentation clock, controlling the precise spatiotemporal patterning of the vertebral column. By counting single transcripts of segmentation clock genes in zebrafish, we show that clock genes have low RNA amplitudes and expression variability is primarily driven by gene extrinsic sources, which is suppressed by Notch signaling. We further show that expression noise surprisingly increases from the posterior progenitor zone to the anterior segmentation and differentiation zone. Our computational model reproduces the spatial noise profile by incorporating spatially increasing time delays in gene expression. Our results, suggesting that expression variability is controlled by the balance of time delays and cell signaling in a vertebrate tissue, will shed light on the accuracy of natural clocks in multi-cellular systems and inspire engineering of robust synthetic oscillators.

Correspondence ertugrul.ozbudak@cchmc.org In Brief Keskin et al. show that segmentation clock transcription levels display low amplitude and high heterogeneity. Clock gene expression variability is primarily driven by gene extrinsic sources, which is suppressed by Notch signaling. Gene expression noise increases along the posteroanterior axis. Spatial gradients of time delays contribute to the noise gradient along the axis.

SUMMARY
Taming cell-to-cell variability in gene expression is critical for precise pattern formation during embryonic development. To investigate the source and buffering mechanism of expression variability, we studied a biological clock, the vertebrate segmentation clock, controlling the precise spatiotemporal patterning of the vertebral column. By counting single transcripts of segmentation clock genes in zebrafish, we show that clock genes have low RNA amplitudes and expression variability is primarily driven by gene extrinsic sources, which is suppressed by Notch signaling. We further show that expression noise surprisingly increases from the posterior progenitor zone to the anterior segmentation and differentiation zone. Our computational model reproduces the spatial noise profile by incorporating spatially increasing time delays in gene expression. Our results, suggesting that expression variability is controlled by the balance of time delays and cell signaling in a vertebrate tissue, will shed light on the accuracy of natural clocks in multi-cellular systems and inspire engineering of robust synthetic oscillators.

INTRODUCTION
Gene expression is inevitably a highly stochastic process due to fluctuations in the complex stoichiometry and reaction kinetics of the biochemical reactions, and it leads to substantial cell-tocell variability (Balá zsi et al., 2011;Elowitz et al., 2002;Kaern et al., 2005;Ozbudak et al., 2002). The resulting phenotypic fluctuations can only be detected and quantified at the single-cell level within isogenic populations. One of the most intriguing questions in science is how developmental pattern formation is executed so robustly despite unavoidable fluctuations in gene expression. This precision necessitates several mechanisms buffering stochastic gene expression. Few studies to date have quantified stochastic gene expression in multi-cellular systems during development (Boettiger and Levine, 2013;Ji et al., 2013;Little et al., 2013;Raj et al., 2010), when buffering the process is critical for the precise and reproducible development of an adult organism, mainly because of technical difficulties posed by quantitative single-cell measurements (Boettiger and Levine, 2013;Ji et al., 2013;Little et al., 2013;Raj et al., 2010).
The study of somitogenesis provides us with an opportunity to investigate the regulation of spatiotemporal precision in pattern formation due to the tight coupling of space and time. The anterior-posterior axis of all vertebrates is patterned as a fixed number of repeating units, the vertebrae. The precursors of vertebrae are derived from somite segments that lay adjacent to the neural tube. Segmentation of somites is dictated by the period of a gene-expression oscillator, called the vertebrate segmentation clock, which is active in unsegmented cells. Oscillatory expression of the Hes/her genes is conserved in vertebrates, and disrupting their oscillatory expression results in vertebral segmentation defects (Pourquié , 2011).
The period of the zebrafish segmentation clock is short: 30 min. At the conclusion of each cycle, a cohort of $200 cells buds from the unsegmented tissue to form a new somite. It remains unknown how segmentation clock factors accumulate to sufficient levels to orchestrate synchronized segmentation of groups of cells. Not only is this process robust and reproducible, but it also undergoes a set number of cycles in a given organism; in zebrafish, segmentation repeats 33 times to form the full-length body axis. Given these robust characteristics, the amplitude of oscillations should be tightly controlled. There is good understanding of the genetic circuitry of the segmentation clock (Pourquié , 2011). However, the field is lacking quantitative measurements of (1) the amplitude of clock gene oscillations, (2) the variability (noise) in clock gene expression between concurrently oscillating cells, (3) how this noise is suppressed in wild-type embryos but unrestrained in certain mutants, and (4) how the magnitude of variability differs spatially among phase-linked cells across a tissue.
To address these questions, we counted RNA molecules transcribed by two master segmentation clock genes (her1 and her7) using single-molecule fluorescence in situ hybridization (smFISH) (STAR Methods) in intact zebrafish embryos. We found low amplitudes and high noise of her1 and her7 transcription in wild-type embryos. Cell-to-cell variability of clock gene expression is dominated by gene extrinsic noise. In Notch signaling mutants, amplitudes of oscillations decreased, while variability increased. Strikingly, transcriptional noise increased spatially from the posterior progenitor zone toward the anterior segmentation zone in both wild-type and Notch signaling mutant embryos. By computational modeling, we showed that the spatially increasing profile of gene expression noise can be recapitulated (B and C) A single z-section of a smFISH picture in a flat-mounted PSM from a wild-type embryo. Cell membrane (green in B), her7 mRNA (red in C) and nuclear staining (blue in B and C). PSM images are rotated 90 degree compared to that in (A). (D and E) The images are zoomed on a highexpression stripe (orange square) (D) or a lowexpression stripe (green square) (E). (F) Tissue is divided in single-cell-wide disks along the axis corresponding to different oscillation phases. Cells containing RNA higher or lower than an arbitrary threshold are plotted as red or gray circles, respectively. Top is left half of PSM, and bottom is right half of PSM. (G) RNA levels from left half of PSM are plotted along the posterior-to-anterior direction. Each dot corresponds to the average RNA number in a spatial phase-binned cell population; error bars indicate 2 SEs. Scale bars mark 30 and 10 mm in (B) and (C) and in (D) and (E), respectively. See also Figure S1 and Table S1. by the spatially increasing gene expression time delays in the clock network.

Segmentation Clock Oscillations Have Low Amplitudes
We combined two strategies in our approach. First, we grouped cells based on their oscillation phases. The oscillation period of segmentation clock genes increases smoothly along the posterior-toanterior (tail-to-head) direction of the presomitic mesoderm (PSM) but remains constant along the left-to-right and dorsal-to-ventral axes (Giudicelli et al., 2007;Gomez et al., 2008). The slowdown of oscillations along the posterior-anterior axis causes a phase delay in cells located in the anterior PSM compared to those located in the posterior PSM. As a result, one sees the different phases of the oscillator cycle mapped out in space along the length of the PSM. Hence, two to three kinematic waves of gene expression can be detected at any moment of an oscillation cycle ( Figures 1A-1C). In other words, all cells located at the same posterior-anterior position in a 2-dimensional, single-cell-wide cross-section are in the same phase of oscillations (Mara et al., 2007;Ozbudak and Lewis, 2008;Riedel-Kruse et al., 2007). Exploiting this unique spatialto-temporal correlation of the segmentation clock, we collected precise, quantifiable data describing the system in detail (Figure 1). Second, we quantified mRNA numbers in single cells in an intact tissue. We used smFISH to quantify the number of RNA molecules in phase-grouped cells along the posterior-anterior axis (Gross-Thebing et al., 2014;Wang et al., 2012). Single cells in the intact tissue were distinguished by jointly analyzing nuclear (DAPI) and cell membrane (membrane-localized GFP) markers ( Figures 1B and S1). Serial sections of fluorescent images of flat-mounted embryos were captured by using a 633 (numerical aperture [NA] = 1.4) objective at 0.240 mm intervals for up to 30 mm (STAR Methods). Thereby, each RNA molecule was detected as a diffraction-limited bright fluorescent dot in the cell ( Figures 1C-1G). The number of clock RNAs per cell did not systematically depend on the z axis, suggesting uniform RNA detection efficiency along the tissue depth ( Figure S1F). To calculate background staining, we used segmented somites, in which expression of clock genes switches off, as a negative control. The background staining was very low (3 and 1 dots per cell for her1 and her7, respectively). To assess the efficiency of RNA detection, we compared the total number of her1 RNA molecules in the PSM of wild-type embryos detected by smFISH to that detected by qRT-PCR (ratio = 1.08 ± 0.37) (STAR Methods). These results validated that a single fluorescent dot in smFISH images corresponds to a single RNA molecule.
We carried out a two-color smFISH experiment to count the number of her1 and her7 RNAs simultaneously in 10-to 14-somite-staged wild-type embryos. In total, we used 18 wild-type embryos. Each embryo captured a particular snapshot of the kinematic waves of the segmentation clock and hence displayed on-off stripes at different positions. The amplitude of oscillations was calculated by subtracting the RNA counts at each spatial position of embryos with low-expressing cells from that of embryos with high-expressing cells ( Figure S2A) and taking the average across spatial positions ( Figure 2A). Although prior computational modeling studies assumed larger amplitude oscillations in the segmentation clock (Ay et al., 2013(Ay et al., , 2014, our data showed that her1 and her7 RNAs have average amplitudes of 41 ± 9 and 49 ± 9, respectively ( Figure 2B). The number of clock transcripts is an order of magnitude lower than those of developmental genes controlling patterning and morphogenesis in Drosophila and housekeeping genes in mammalian cells (Boettiger and Levine, 2013;Little et al., 2013;Padovan-Merhar et al., 2015). It is remarkable that embryos can robustly accomplish segmentation by transcribing two master oscillating genes at such low levels, where they would be more vulnerable to gene expression noise (Balá zsi et al., 2011;Elowitz et al., 2002;Ozbudak et al., 2002). Next, we investigated the magnitude of expression noise in the clock genes.
Gene Extrinsic Noise Dominates the Total Cell-to-Cell Expression Variability The two master segmentation clock genes (her1 and her7) are adjacent but separated by a common regulatory sequence and transcribed in opposite directions ( Figure 2C). These genes are paralogous duplicates and code for transcriptional repressors that function as hetero-or homo-dimers. Expression of clock genes is repressed by their protein products, which create a negative feedback loop that is thought to be the pacemaker of the segmentation clock (Harima et al., 2013;Lewis, 2003). her1 and her7 have similar transcriptional time delays (Hanisch et al., 2013) and RNA half-lives (Giudicelli et al., 2007). Transcription of her1 and her7 is therefore concomitant (Gajewski et al., 2003;Oates and Ho, 2002), and our preceding results showed that they have comparable amplitudes ( Figure 2B). Intrinsic fluctuations are those due to the randomness inherent to transcription; being random, they should affect the transcription of her1 and her7 independently, producing uncorrelated variations in respective mRNA levels ( Figure 2D) (Elowitz et al., 2002). Other molecular species in the cell, e.g., RNA polymerases and upstream transcriptional regulators, are gene products and therefore will also vary over time and from cell to cell. These variations cause correlated fluctuations in the expression of clock genes and are defined as extrinsic noise (Elowitz et al., 2002). To separate the contributions of intrinsic and extrinsic sources in the expression variability of these two clock genes, we took an approach similar to earlier works ( Figure S3; Supplemental Experimental Procedures) (Elowitz et al., 2002;Rhee et al., 2014).
We first calculated expression variability among phase-grouped cells (single-cell-diameter spatial slices) (Figures 2E and 2F) and then further grouped the slices into five bins based on mean clock mRNA levels ( Figure 2G). Our results showed that the total expression variability of segmentation clock genes, as described by CV 2 ([SD/mean] 2 ), is 0.49 (i.e., CV = 0.7) at the lowest transcription state and 0.20 (i.e., CV = 0.45) at the highest transcription state ( Figure 2G), which is substantially higher than the measurement error ( Figure 2G; Supplemental Experimental Procedures). The total clock gene expression noise ( Figure 2G) is as high as the levels observed in single-cell systems (Golding et al., 2005;Raj et al., 2006;Taniguchi et al., 2010) but much higher than those of developmental genes in Drosophila (Boettiger and Levine, 2013;Little et al., 2013), likely due to lower transcript levels ( Figure 2B). The gene extrinsic noise contributes most to total variability in transcription (p < 0.001), while the gene intrinsic noise is only high at low transcript levels and decreases proportionally to the total clock transcript level (Figure 2G). Slight differences in the promoter strengths of her1/ her7 genes and transcriptional time delays (Hanisch et al., 2013) and half-lives of her1/her7 mRNAs (Giudicelli et al., 2007) might have slightly increased gene intrinsic noise in our calculations. However, our conclusion that extrinsic noise dominates clock transcriptional variability holds true. In yeast, noise in gene expression is primarily extrinsic in origin (Becskei et al., (F) Total noise (CV 2 = [SD/mean] 2 ) versus mean levels of total her (her1 + her7) RNA is plotted. Each blue dot represents the values in a single slice. All slice data from all wild-type embryos are pooled. (G) Total noise split into intrinsic and extrinsic components for each slice and all slice data reported in (F) grouped into five bins based on mean her RNA levels. The intrinsic noise decreases with the average RNA levels (magenta). Extrinsic noise (green) is the dominant component of total expression noise (blue), p < 0.001 (low and high expression). Total measurement error due to RNA counting, cell segmentation, and phase grouping of cells is plotted as a baseline (dashed black curve). y axis is noise; x axis is mean levels of total her (her1 + her7) RNA in grouped slices. The graph is in log-log scale; error bars indicate 2 SEs. See also Figures S2 and S3. (B) Total RNA amplitudes decreased in deltaC À/À (green) and deltaD À/À (red) mutants compared to wild-type embryos (purple) and in DAPT-treated embryos (blue) compared to DMSO-treated control embryos (pink) (p < 0.001).
(legend continued on next page) 2005; Colman- Lerner et al., 2005;Raser and O'Shea, 2004;Volfson et al., 2006); however, the intrinsic noise is found to dominate for hb patterning gene in Drosophila (Little et al., 2013). Cell cycle is less likely to contribute to extrinsic noise, because the segmentation clock period is much smaller than cell cycle and only a small population of cells experience mitosis during an oscillation cycle in the zebrafish PSM (Horikawa et al., 2006). We next assessed whether variation in cell size could underlie the dominant extrinsic noise. We reanalyzed our data by normalizing transcript counts by cell volume. The results showed that cell volume has only a minor contribution on the expression noise and extrinsic noise is the dominant portion of clock variability ( Figure S2B). Because her1 and her7 genes are located adjacent to each other, common upstream factors (Becskei et al., 2005;Volfson et al., 2006) or local chromatin modifications are likely the main factors resulting in dominant extrinsic noise, as shown in yeast (Becskei et al., 2005), but one cannot exclude variations in epigenetic and metabolic states.

Notch Signaling Buffers Mostly Gene Extrinsic Noise
After finding that extrinsic noise dominates the gene expression variability in the segmentation clock, we investigated the mechanism that might limit its magnitude. In Drosophila, noise in nascent transcription of hb gene is as high as that of the zebrafish segmentation clock genes, but noise in Drosophila patterning genes is averaged out by spatial averaging in syncytial embryo and by temporal averaging due to longer half-life of RNA (Little et al., 2013). However, spatial averaging is not possible in the vertebrate PSM cells, because they are not syncytial, and the short half-lives of RNAs (3-8 min) (Giudicelli et al., 2007) prevent simple temporal averaging.
Oscillatory membrane-bound DeltaC ligands activate Notch receptors in neighboring cells; hence, Notch signaling acts at a short distance ( Figure 3A). The activated intracellular domain of the Notch receptor, in turn, activates transcription of her-family genes. This short-distance coupling synchronizes oscillation phases among cells located in the same anterior-posterior position in the PSM (Delaune et al., 2012;Horikawa et al., 2006;Jiang et al., 2000;Mara et al., 2007;Ozbudak and Lewis, 2008;Riedel-Kruse et al., 2007). Based on experimental observations and computational modeling, Notch signaling has been hypothesized to minimize gene expression noise (Delaune et al., 2012;Horikawa et al., 2006;Jiang et al., 2000;Mara et al., 2007;Ozbudak and Lewis, 2008;Riedel-Kruse et al., 2007), but gene expression noise has not been quantified so far. To quantify noise, binary descriptions such as ''synchronous versus asynchronous (salt-and-pepper)'' are insufficient; one must count single molecules (Raj and van Oudenaarden, 2009) and quantify their distributions in cells that are grouped according to their oscillation phases. To investigate the impact of Notch signaling-mediated cell-to-cell coupling on expression noise, we carried out two-color smFISH experiments in two mutant lines-deltaC À/À (J€ ulich et al., 2005) and deltaD À/À (Holley et al., 2000)-in which Notch signaling is disrupted and stripy expression patterns are lost due to desynchronized oscillations. The amplitude of oscillations and mean RNA levels decreased (ranging from 23% to 54%) in two mutants compared to wildtype embryos (p < 0.006) ( Figures 3B and S3). Because transcriptional noise depends on mean RNA levels ( Figure 2G), we compared noise between mutants and wild-type embryos at similar mean RNA levels. The mutants lacked higher expression groups due to reduced clock transcription. Total expression noise increased (more than 30%) in two mutants compared to wild-type embryos due to extrinsic (p < 0.001), but not intrinsic, noise ( Figures 3C and 3D). Extrinsic noise was elevated predominantly in the lower expression states. However, the mutant data reflect the steady-state response and might miss the initial outcome of loss of Notch signaling. To investigate this issue further, we transiently blocked Notch signaling by treating embryos with the g-secretase inhibitor N-[N-(3,5-difluorophenacetyl)-L-alanyl]-S-phenylgly cine t-butyl ester (DAPT) or DMSO (as a control) for 1.5 hr. Earlier work showed that long-term DAPT treatment recapitulated segmentation defects in mutants (Mara et al., 2007;Ozbudak and Lewis, 2008;Riedel-Kruse et al., 2007). Two-color smFISH experiments showed that blocking of Notch signaling for only 1.5 hr did not disrupt the stripe expression patterns ( Figures 3E and 3F) as reported previously, suggesting that the global phase synchrony is not disrupted with short-term DAPT treatment (Mara et al., 2007;Ozbudak and Lewis, 2008;Riedel-Kruse et al., 2007). However, the amplitude of oscillations decreased (28%, Figures 3B and S3) while the expression noise slightly increased in the DAPT-treated embryos compared to DMSO-treated embryos. As expected, the elevation of noise in DAPT-versus DMSO-treated embryos (Figures 3G and 3H) is much lower than that in mutants versus wildtype embryos (Figures 3C and 3D) due to the transient nature of drug treatment. Increased noise would then accumulate in time to cause both amplitude variability and phase drifting among neighboring cells and result in higher variability in mutant (deltaC À/À and deltaD À/À ) embryos. These results affirmed the long-standing hypothesis that the cause of segmentation defects in Notch signaling mutants is untamed gene expression noise in the segmentation clock genes (Jiang et al., 2000) and discovered that the main function of Notch signaling is to dampen extrinsic noise.
Although Notch signaling is an activator of clock transcription ( Figure 3B), total noise did not primarily increase due to the decreased clock transcription in Notch mutants for the following (C and D) Extrinsic noise (C) increased in deltaC À/À (green) and deltaD À/À (red) mutants compared to wild-type embryos (purple) (p < 0.001 in both lower and higher groups). Intrinsic noise levels (D) are shown in deltaC À/À (green) and deltaD À/À (red) mutants and wild-type embryos (purple). x axis is mean levels of total her (her1 + her7) RNA in grouped slices. Clock transcription decreased in notch mutants. Therefore, slices could only be grouped for the lower expression bins compared to wild-type data. (E and F) A single z-section of a smFISH picture for her7 mRNA (red) in DMSO-treated embryos (as a control, E) and DAPT-treated embryos (F). The stripy expression pattern of her7 mRNA is not disrupted after 1.5 hr DAPT treatment (F). Scale bars mark 30 mm. (G and H) Extrinsic (G) and intrinsic (H) noise levels in DAPT-treated embryos (blue) and DMSO-treated control embryos (pink). Error bars indicate 2 SEs. See also Figure S3.
reasons. First, gene intrinsic noise displays an inverse relationship with mean expression levels (Elowitz et al., 2002;Kaern et al., 2005) ( Figure 2G). However, total noise increased in mutants due to increased extrinsic, but not intrinsic, noise ( Figures  3C and 3D). Second, the amplitude of oscillations reduced twice in deltaD À/À mutant compared to deltaC À/À mutant, but intrinsic noise was not higher in deltaD À/À mutant than in deltaC À/À mutant. Altogether, these results show that total noise is elevated in Notch pathway mutants due to both desynchronized oscillations and decreased amplitudes, and surprisingly, these two factors increase only extrinsic, not intrinsic, noise. Extrinsic noise is the dominant component of the transcriptional variability in wild-type embryos ( Figure 2G), which implies that transcriptions of her1 and her7 genes fire nearly simultaneously in each PSM cell. When Notch signaling is disrupted, extrinsic noise is elevated further (Figures 3C and 3D). This result shows that Notch signaling is not the source of co-transcription of two clock genes but rather actively tames gene extrinsic noise. We hypothesize that high extrinsic noise is likely due to the cell-to-cell variability in either the levels of other transcriptional regulators or the epigenetic state of the chromosomal locus.
Gene Expression Noise Increases from the Progenitor to the Segmentation and Differentiation Zone Afterward, we assessed how noise in clock gene expression varies in space from the posterior PSM (progenitor zone) to the anterior PSM (where cells initiate segmentation) by pooling the noise values at each spatial position from all wild-type embryos ( Figure 4A) (see STAR Methods for details). Strikingly, we found that expression noise increases toward the anterior PSM in wildtype embryos (p < 0.001) ( Figure 4B). In Notch signaling mutants (deltaC À/À and deltaD À/À ), the spatial profile of gene expression noise does not flatten but rather elevated throughout the PSM (p < 0.001) ( Figure 4B). Because total noise (CV 2 ) depends on average expression levels ( Figure 2G), we assessed the spatial change of noise by grouping data based on mean expression levels (see STAR Methods for details). Transcriptional noise is higher in cells located in anterior PSM than those in posterior PSM at comparable mean RNA levels (p < 0.001) (Figures 4C, 4D, and S4; STAR Methods). The results showed that the total noise is minimized by Notch signaling throughout the PSM in wild-type embryos but is boosted in Notch signaling mutants, and Notch signaling does not play a role in the anteriorly increasing trend of noise in wild-type embryos (Figure 4).
These results could be perceived as unexpected, because cell-to-cell variability should be minimized for cohorts of cells to segment robustly at the end of anterior PSM. This paradox can be explained by our earlier results, in which we have shown that the segmentation clock functions primarily in cells located in the posterior PSM. The segmentation clock relays its periodic information to downstream genes in the middle of the PSM (Giudicelli et al., 2007); thereafter, the clock is not needed for the segmentation process to be completed, although it continues to be expressed in the anterior PSM. It has been puzzling why expression of the segmentation clock became functionally irrelevant for cells located in the anterior PSM. Our results demonstrating increased gene expression variability of clock genes in anterior PSM reveal why cells could not rely on noisy clock expression in the anterior PSM. Collectively, our results show that there is selection pressure to minimize noise in clock gene expression only in the posterior PSM ( Figure 4B). These results are the first demonstration of spatial variation of gene expression noise during vertebrate development.
We further investigated the mechanism driving spatial gradient of clock genes' expression noise. A candidate mechanism is the spatially increasing effective gene expression time delays along the posterior-anterior direction (Ay et al., 2014). We built a simple stochastic computational model to test the impact of time-delay gradients on the expression noise in the tissue (Figures 4E and S5;Supplemental Experimental Procedures). In our model, we simulated synchronized oscillations in two neighboring cells. Cell-autonomous oscillations are generated due to an intracellular negative feedback loop (with an effective time delay of t x ), and the clocks are synchronized due to an intercellular positive feedback loop (with an effective time delay of t y ). We then repeated simulations by increasing the time delays gradually up to 4.5-fold-a value that matches to the experimentally quantified increase in effective time delays along the PSM (Ay et al., 2014). The simulations recapitulated the same spatially increasing trend in gene expression noise as time delays are increased along the tissue ( Figure 4F).

DISCUSSION
Many studies of stochastic gene expression have been performed with synthetic regulatory networks or promoter reporters for natural networks by using long-lived GFP mRNA and protein. However, many mRNAs coding for critical transcription factors governing embryonic development and adult homeostasis are short lived (Schwanhä usser et al., 2011). Here, by studying short-lived segmentation clock RNAs (t 1/2 = 3-5 min) (Giudicelli et al., 2007), we investigated gene expression noise for the first time at the fastest dynamic scale (30-min period) in an intact vertebrate tissue.
Earlier work evaluated her/hes-family gene expression in embryos by simple binary scoring (signal present versus absent) (Delaune et al., 2012;Horikawa et al., 2006;Jenkins et al., 2015;Jiang et al., 2000;Mara et al., 2007;Ozbudak and Lewis, 2008;Riedel-Kruse et al., 2007), by arbitrary fluorescence units instead of absolute molecular counting (Webb et al., 2016), or in cell culture in which spatial information and Notch coupling were lost (Phillips et al., 2016;Webb et al., 2016). Here we used singlemolecule counting in an intact vertebrate tissue to analyze the amplitude and variability of the segmentation clock genes in different genetic backgrounds and spatial positions. Our results showed that the amplitudes of transcriptional oscillations during somitogenesis are very low compared to RNA levels of developmental genes in other systems (Boettiger and Levine, 2013;Little et al., 2013) and extrinsic noise contributes the most to the high variability of clock gene expression (Figure 2). The period of oscillations is short ($30 min), and transcriptional time delays make up a significant portion of the period (Giudicelli et al., 2007). This would restrict the number of proteins that can be translated from a transcript in a given clock cycle. Unlike constitutively expressed and stable proteins, the segmentation clock proteins are short lived (Ay et al., 2013). Therefore, clock protein levels closely follow mRNA levels. Furthermore, clock RNAs and proteins establish a transcriptional negative feedback loop. All these features make clock mRNA levels the best reporters of clock protein levels. Here, we hypothesize that a major portion of tran-scriptional extrinsic noise is due to variability in clock protein levels. Although clock proteins have not been counted yet, published data revealed extensive cell-to-cell variability in Her1-Venus reporter levels (Delaune et al., 2012). In this study, we (B) The total her expression noise increased in deltaC À/À (green) and deltaD À/À (red) mutants compared to that in wild-type (purple) throughout the PSM, p < 0.001 (wild-type and deltaC À/À in posterior and anterior), p < 0.001 (wild-type and deltaD À/À in posterior and anterior). x axis is spatial positions. (C) Representative sketches of n embryos at different oscillation snapshots are shown. Noise values from posterior and anterior PSM of all embryos are grouped in the first and second groups, respectively. (D) Expression noise is higher in cells located in the posterior versus the anterior PSM, p < 0.001. Expression noise (y axis) decreases as expression levels (x axis) increases, as shown in Figure 2G. However, at a given expression level, cells located in the anterior PSM have higher expression noise for the clock genes than cells located in the posterior PSM. The spatially increasing trend of expression noise (B) is independent of mean expression levels. have investigated the precision in the amplitude of oscillations. Investigating the precision in the timing of oscillations (Morelli and J€ ulicher, 2007) would require carrying out single-molecule RNA or protein counting in real time. The levels of Her1/Her7 proteins should be above a threshold to repress many target genes, including deltaC; precision in their levels is therefore important for the successful segmentation of somites.
The timely and precise progression of a genetic program is critical to achieving reproducible embryonic development. Unlike other systems (Little et al., 2013;Shah and Tyagi, 2013), in somitogenesis, gene expression variability is not averaged in time or space. Here, we showed that vertebrate embryos use Notch signaling-mediated cell-to-cell coupling to minimize expression variability throughout the PSM (Figure 3). In the absence of Notch signaling, increased gene expression noise ( Figure 3) desynchronizes clock oscillations and results in segmentation defects (Mara et al., 2007;Ozbudak and Lewis, 2008;Riedel-Kruse et al., 2007) and congenital scoliosis in patients (Pourquié , 2011). Lastly, we found that transcription noise increases spatially as cells approach the segmentation zone, likely due to spatially increasing gene expression time delays but independent of Notch signaling (Figure 4). Because time delays are inevitable in the expression of every gene, understanding the role of time delays in expression noise is significant. We showed that the spatial gradient of effective time delays in the PSM (Ay et al., 2014) likely contributes to the spatially increasing expression noise of the segmentation clock genes along the axis ( Figure 4F).
Oscillations are widespread in biological systems. Notch signaling plays critical roles in controlling the switch from proliferation to differentiation in almost every tissue throughout the metazoan. Hes/her family genes are direct targets of Notch signaling. The segmentation clock was the first discovered developmental oscillator (Palmeirim et al., 1997). However, Hes/Her protein levels also oscillate in neural progenitors, embryonic stem cells, and ovarian cells, thereby controlling the switch from proliferation to differentiation in several tissues (Kobayashi and Kageyama, 2014). Furthermore, Hes/Her proteins are highly expressed in several types of human tumors, while their inhibition has been shown to restore differentiation (Kobayashi and Kageyama, 2014;Liu et al., 2015;Sang et al., 2010). We anticipate that future studies may lead to developing new ways to control stem cell proliferation and differentiation in various tissues, developing therapies against certain cancer types, understanding the precision of other natural oscillators, and engineering synthetic oscillators.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

DECLARATION OF INTERESTS
The authors declare no competing interests.  (Cooper et al., 2005) transgenic line was used as wild-type, and dlc tw212b/tw212b (J€ ulich et al., 2005) and dld tr233/tr233 (Holley et al., 2000) were used as Notch signaling mutants. Temporal loss of function of Notch signaling was accomplished by treating embryos with 100 mM of the g-secretase inhibitor, N- [N-(3,5-difluorophenacetyl)-L-alanyl]-S-phenylgly cine t-butyl ester (DAPT) while DMSO is used as a control (Ozbudak and Lewis, 2008). Fish were bred and maintained at 28.5 C on a 14-10 hr light/dark cycle.

METHOD DETAILS smFISH and imaging
We have used the protocol developed by Advanced Cell Diagnostics (Wang et al., 2012), which allows for simultaneous detection of two differentially labeled transcripts in zebrafish (Gross-Thebing et al., 2014). 0.9 nL of membrane-localized GFP RNA was injected in one cell stage mutant embryos. Afterward, embryos were incubated at 23 C until they reach 10-14 somite stage, then embryos are fixed in 4% PFA in PBS for 1 hr at room temperature (RT), washed with 0.1% PBSTw (0.1% Tween-20 in PBS), then with gradually increasing methanol (MeOH) concentration (50% MeOH -50 PBSTw, 100% MeOH), each wash for 5 min. Dechorionation was performed when they were in 50% MeOH -50% PBSTw solution. Afterward, embryos were stored in 100% MeOH at À30 C overnight (O/N). The next day, embryos were air-dried for 30 min at RT and processed with Pretreat 3 (Advanced Cell Diagnostic RNAScope Pretreatment Reagents, 320842) for protease digestion for 20 min at RT. Then, they were washed with 0.01% PBSTw (0.01% Tween-20 in PBS) for 3 times each for 5 min. Probe-Dr-her7 and Dr-her1-LE2-C3 probes were mixed in 50:1 ratio, respectively, and probe mix was warmed at 40 C for 10 min in the oven, cooled down to RT, and added on embryos. Embryos were incubated in hybridization oven at 40 C for O/N. The next day all washing and signal amplification steps with RNAscope Fluorescent Multiplex Detection Reagents (Advanced Cell Diagnostics 320851) were performed according the RNAscope Protocol for Zebrafish (Gross-Thebing et al., 2014), with only few changes. Two additional washing steps were added both after the removal of probe and after the pre-amplifier hybridization step to decrease the background signal (5 times for 15 min) and each reagent was used as 1 drop. For labeling Amp4 Alt-B reagent was used. After the labeling step with Amp4, embryos were fixed in 4% PFA in PBS at 4 C for O/N. The following day, embryos were washed with 1% PBSTX (1% Triton X-100 in PBS) for 3 times each for 5 min, and permeabilized in 1.5% Triton X-100 for 1 hr at RT. Embryos were then incubated in blocking buffer (1% Triton X-100, 2%BSA, 5% Goat Serum) for 2 hr at RT. Then, embryos were incubated with Rabbit anti-GFP primary antibody (Life Technologies A6455) in blocking buffer (1:100) at 4 C for O/N. The following day, embryos were washed with 1% PBSTX 3 times for 10 min followed by a wash with blocking buffer for 10 min. Then, embryos were incubated in blocking buffer with 1:400 Hoechst trihydrochloride, trihydrate (Invitrogen 33342) and 1:200 Goat anti-Rabbit IgG Alexa 488 Secondary Antibody (A-11034, Life Technologies) at 4 C in dark for O/N. Then, embryos were washed with 0.2X SSCT (0.01% Tween-20 in 0.2X SSC) for 15 min at 4 C and fixed in 4% PFA in PBS. Embryos were mounted in 0.2X SSCT solution. ProLong Gold antifade reagent (Life Technologies P36934) was used to prepare slides. Images were captured by using a 63X (NA = 1.4) objective at the Zeiss Imager Z2. Serial sections of fluorescent images were taken at 0.240 mm intervals for up to 30 mm with AxioCamMRm camera, apotome and Axiovision software 4.8.2 Release. Images of single embryos were tiled along the x-y axis and stitched by the Axiovision software.
Quantifying RNA numbers in single cells Imaris 8.1.2 Cell Module was used to analyze the images. First, the tissues surrounding the PSM (notochord, neural tube and lateral plate mesoderm were filtered by using the surface tool. Manual surface creation tool was selected to mark the regions covering undesired tissues and the voxels inside were set to zero. Individual PSM cells were identified by using the cell tool. Nuclei diameter threshold was set to 3 mm with smoothing width of 0.3 mm, background subtraction sphere diameter of 1.2 mm and split nuclei by seed points. Then, number of voxels filter was used to select the seed points and an intensity threshold was used for background subtraction, finally nucleus number of voxels was set to select single nucleus. Afterward, detect cell boundary from cell membrane option with 0.25 mm membrane width was selected. The cell volume filter was set between 150-450 mm 3 . Finally, spots tool was used to count the total number of RNA molecules in the PSM. Estimated diameter was set to 0.5 mm and background subtraction option was selected. Quality score filter was set to separate RNA dots from background. The cell-membrane signal (GFP) got weaker in the deepest z-layers, and thus we could not successfully separate cells and collect data on the ventral-most part ( Figure S1). Therefore, we have only included data from cells that we could successfully segment by image analysis. The number of clock RNAs per successfully segmented cell do not systematically depend on z axis ( Figure S1F).

qRT-PCR
Full-length her1 mRNA standard was synthesized from her1 pCS2+ clone. The plasmid was linearized by Not1 and SP6 mMessage mMachine (Life Technologies AM1340) was used for in vitro transcription. RNA was purified with Zymo Research Quick RNA MiniPrep Kit (R1054). The loss of RNA during column purification was calculated by measuring the concentration of RNA by Qubit (Life Technologies) both before and after the column purification. Tails of five 13-14 somite-staged embryos were dissected and pooled in L15 media under a dissecting microscope and purified with Zymo Research Quick RNA MiniPrep Kit. Five biological replicates of this experiment were carried out. RT reactions were performed by using SuperScript IV Reverse Transcriptase (Invitrogen 18090010), RNaseOUT (Invitrogen 10777-019), dNTP (Roche 05 081 955 001) and the reverse primer: (5 0 -ACG TCT CGA GTC ACC AGG GTC TCC ACA AAG GCTG-3 0 ). q-PCR reactions were performed by using the forward primer: 5 0 -TGG AAG AAC TGC GAA CGC TT-3 0 , reverse primer: 5 0 -CGGAGGTTTTGGATCATGCG-3 0 and SYBR Select Master Mix (Life Technologies 4472908) with the Applied Biosystems StepOnePlus Real-Time PCR System. Standard RNA dilution series were prepared by 10-fold dilutions and the efficiency of primers was calculated by fitting a standard curve in the StepOne Software v2.3 (Slope: À3.32, R 2 : 0.98, Y-intercept: 33.423, Efficiency: 100%). The numbers of total her1 mRNA in tail samples were calculated by comparing their C T values to that of the standard her1 RNA sample. We have detected 82264 ± 10816 her1 mRNAs per single embryo.
Calculating the efficiency of RNA detection by smFISH We counted the total RNA numbers in the whole PSM by smFISH and compared the results to that obtained by qRT-PCR. Spots tool was used to count the total number of her1 RNA molecules in the PSM of 18 wild-type embryos. Estimated diameter was set to 0.5 mm and background subtraction option was selected. Quality score filter was set to above 70. We subtracted the nonspecific smFISH staining as follows: Surface tool was used to count the total number of nuclei (cells) in the PSM by using the DAPI channel data.
Smooth tool was selected. Surface area detail level was set to 0.3 mm. Background subtraction with local contrast was used and the diameter of largest sphere which fits into the object was set to 1.2 mm. The threshold for background subtraction was set to 100. Split touching objects was enabled with seed point diameter as 3 mm. The number of voxels filter was used to separate the nuclei from smaller objects. The total number of cells was multiplied with the background smFISH staining and the resultant values were subtracted from that obtained with the spot detection tool to obtain the total number of her1 mRNA molecules in a single embryo: 75862 ± 24061. The ratio of the number of her1 mRNA molecules detected by qRT-PCR versus smFISH is: 1.08 ± 0.37.
Calculating the measurement error in RNA counting by smFISH Control smFISH experiments were carried out in wild-type embryos by using two colors of probes (Dr-her1-LE1-C1 and Dr-her1-LE2-C3) directed against the her1 mRNA. The probe sets are designed in an alternating tiled manner complementary to the her1 mRNA. In total, we have used 8 wild-type embryos and imaged 21-35 sections in each embryo. The number of RNAs in each cell is measured as described in the above sections. Since binding of one probe set to RNA would be independent of the other probe set, the measurement error in counting single RNA molecules can be calculated by using an equation mimicking intrinsic noise: ð1=2Þ<ððher1C1=<her1C1 > Þ À ðher1C3=<her1C3 > ÞÞ 2 > . The measurement error in RNA counting is: 0.03.

Calculating the measurement error in cell segmentation
One of the sources of measurement error would be incorrectly segmenting cells from each other, i.e., failing to assign image voxels unambiguously to a single cell. To estimate the cell segmentation error, we in total selected 38 image sections from 7 different wildtype embryos. We manually segmented 255 cells from this dataset by ImageJ. For each cell, we compared the manually segmented areas in ImageJ with the automatically segmented areas in Imaris. The measurement error of misassigning mRNAs to single cells is calculated as follows: Let x i be an independent and identically distributed (iid) random variable denoting the actual mRNA count in a cell with mean hxi and coefficient of variation squared CV 2 . Assuming each cell has a normalized area of 1, the measurement is given by where a i is an iid random variable denoting the fractional area ignored, x j is the mRNA count in the neighboring cell, and c i is an iid random variable denoting the fractional area gained of the neighboring cell. The mean and variance of a i is given by hai and s 2 a , respectively. Similarly, the mean and variance of c i is defined as hci and s 2 c . Then, the error in the measurement would be Taking the variance of (2) and dividing by hxi 2 quantifies the error in terms of the coefficient of variation squared: We calculated the segmentation error at different total her mRNA levels ( Figure 2G). The measurement error values in cell segmentation ranged between 0.02 and 0.03.

Grouping cells in space based on oscillation phases
The position of each cell and the number of her1 and her7 RNA molecules in each cell were measured in each embryo. Each embryo was divided into left and right halves by visual inspection ( Figure 1F). Within each half PSM tissue, cells were grouped (sliced) based on their oscillation phases, which vary smoothly along the axis (Giudicelli et al., 2007). The slice width was set to 8 mm, which corresponds to the diameter of cells in the PSM. The angles of slices were fixed to the angle of expression stripes of the segmentation clock in the PSM. The angle of expression stripes of the segmentation clock gene changes incrementally along the posterioanterior direction in the PSM. We first measured the stripe angles along the axis in all wild-type and DAPT-/DMSO-treated samples. We then fitted an equation to the data ( Figure S1G). We obtained the following equations: wild type angle = 0:039 Ã Distance + 44:23, DMSO angle = 0:057 Distance + 43:45, and DAPT angle = 0:123 Distance + 21:43, where distance (mm) is measured from the tail end of embryos. Later on, we used these equations to incrementally change the angles of slices along the PSM in all embryos in each background. For deltaC and deltaD mutant embryos, we have used the angle function of the wild-type embryos. We assigned each cell into a slice when the center of the cell is located within a spatial slice.
Calculating the measurement error in phase grouping of cells The calculation of expression noise depends on assignment of cells to correct oscillation phase bins. Part of the extrinsic noise could be due to slight errors in the assessment of cell positions rather than coregulation. To address this issue, the anterior stripe expression angles for wild-type embryos are calculated manually by two different experimentalists. We calculated two total expression noise values at different total her mRNA levels by using these two different angle sets. The measurement error is set to be the difference between total expression noise values obtained by using two independent measurements of expression angles. The measurement error in phase grouping of cells ranged between 0.01 and 0.05.

Calculating the total measurement error
We added the measurement errors due to RNA counting, cell segmentation and phase grouping of cells to obtain the total measurement error in our experimental and analysis pipeline at five total her mRNA levels ( Figure 2G). The total measurement errors were between 0.066 and 0.098. We plotted the total measurement error in Figure 2G as a baseline to biological variability (expression noise) that can be unambiguously assayed.

QUANTIFICATION AND STATISTICAL ANALYSIS Heatmaps
For visual aid, we created heatmaps for each embryo by grouping the cells into two groups: cells with low or high mRNA numbers ( Figure 1F). The threshold for this grouping was chosen as follows. First, the minimum (N) and maximum (M) mRNA levels per cell were determined in the embryo. Then the threshold was chosen as N+0.3*(M-N) by visual inspection. The cells below this threshold were assigned as low expression cells, and the rest were assigned as high expression cells. The heatmaps were created by plotting high expression cells with bold colors and low expression cells with light colors. Slice boundaries were shown in each heatmap (as in Figure 1F). The details of the data analysis and computational modeling can be found in Supplemental Experimental Procedures.

Subtracting nonspecific staining
In each genetic background, we quantified the nonspecific smFISH staining in the segmented somites, where her1 and her7 genes are not expressed. Mean and variance of nonspecific staining were given in Table S1.
Mean nonspecific staining levels were subtracted from each cell's mean her1 and her7 mRNA levels. Afterwards, cells with negative her1 or her7 mRNA levels were removed from the data. Slices with fewer than three cells after background subtraction were also removed from the data.

Spatial gene expression dynamics
Raw data was processed differently for this analysis. If a slice had negative mean of either her1 or her7 mRNA after background subtraction, the negative mean expression level was set to zero and the slice was not eliminated. For left and right halves of each embryo, her1 and her7 mean mRNA levels (y-axis) at each spatial location (slice) were plotted (x-axis). Spatial location 0 corresponded to the posterior (tail) end of embryos ( Figure 1G).

Amplitude measurement
Slices with negative mean mRNA level were not eliminated but rather mRNA levels were set to zero. Amplitudes of her1 and her7 gene expression levels were calculated as follows: Mean mRNA levels of each slice from all embryos were used to calculate amplitude for each genetic background. In this approach, we first combined the data from individual embryos in each genetic background as follows. First, we aligned the posterior ends of all embryos. Then, we combined the data from the ith slices from both left and right halves of all embryos. This iterative process was performed until one of the embryos ran out of slices. By this way, we collected 36, 28, 36, 18 and 24 data points for each cell position from 18 wild-type, 14 deltaC -/mutant, 18 deltaD -/mutant, 9 DAPT-treated, and 12 DMSO-treated embryos, respectively. To reduce the effect of outliers: a) We grouped data from every five consecutive spatial locations.
Within each spatial group, mRNA levels of all slices were ranked from largest to smallest. b) Mean of the top 10% of the group was assigned as the peak and mean of the bottom 10% was assigned as the trough. The difference between the peak and trough in each group was defined as the amplitude at that location. We tested the effect of binning on amplitude values by varying the number of spatial groupings from 4 to 8 and top/bottom groupings from 5% to 10%. These modifications did not change the amplitude values more than 16%. The spatial profile of total her mRNA amplitude is plotted in Figure S2A. We then averaged the amplitude values across the space and plotted the average amplitudes in Figure 2B.

Noise formulas
The intrinsic and extrinsic noise levels for each slice were computed using the equations below.
Total noise was calculated as the sum of intrinsic and extrinsic noise.

The contribution of differences in cell volumes on expression noise
We calculated the concentration of mRNAs in each cell by dividing the number of mRNAs in each cell by the cell volume. Later, we normalized the resulting values to report the number of mRNAs per 4 pl volume. After this normalization, the noise plots reported in Figure   2G and Figure S2B had comparable ranges on the x-axis.

Spatial variation in noise
The CV 2 data from all slices in all embryos in each genetic background was combined as described in the "Amplitude measurement" section. CV 2 values of the ith slices from both left and right halves of all embryos were combined. The mean and standard error of CV 2 values were plotted at different spatial positions ( Figure 4B).
To control the dependency of spatial CV 2 values on expression levels, we carried out two further analyses: 1) The data were grouped into three categories depending on its expression level: low, medium and high. Slices were first ranked from the highest to the lowest gene expression levels.
Then, they were divided into three equally sized groups: bottom 33.33% (low expression), middle 33.33% and top 33.33% (high expression). For each spatial location, three different CV 2 values (at low, medium and high mRNA levels) and their standard errors were calculated by using the sum of her1 and her7 mRNA levels at each genetic background (wild-type, deltaC -/and deltaD -/-; Figure S4A).
2) The embryos' slices are placed into 8 groups based on their mean total her expression levels. Slices with expression levels from 0 to 15 formed the first group, and slices from 105 to 120 formed the last group. The slices above 120 are considered to be outliers, and eliminated from the analysis. Afterwards, for each expression group, the data coming from slices located in the posterior 33% and anterior 33% of the PSM are selected. Then, the mean plus two standard error of expression level, and mean plus two standard error of CV 2 in posterior and anterior slices are calculated. The figure summarizing this analysis is provided in Figure 4C, D.

Statistical analysis
One-way ANOVA was conducted to determine the statistical significance of the differences between amplitude and noise in different genetic backgrounds and/or spatial locations. Normality was assessed by visual inspection of histograms and normal Q-Q plots.
Homogeneity of variances was assessed by comparing standard deviations of different groups.

Computational Model
We built a computational model to simulate synchronized oscillation in two neighboring cells. In this model, a cell (Cell a) produces a protein > from a constitutive gene as shown in Figure S5A. This protein activates the production of two other proteins > and > . The protein > inhibits the production of the protein > after a time-delay B , while the protein > activates the production of in the neighboring cell (Cell b) after another delay C . The same analogy is true for the neighboring cell. Overall, the production rate of protein > at time is given by where > ( G ) denotes the count level of > ( G ) at time − B ( − C ). Further E , B , and C are basal production rate, negative feedback strength, and positive feedback strength, respectively. Finally, protein decays with rate .
The biological representation of the aforementioned time-delays is as follows: The timedelays in Figure S5A are implemented by considering that the proteins and will be activated after some time. Hence, we consider that the protein > produces two intermediate molecules K Assuming that the timing of first-order reactions are similar and exponentially distributed, the conversion process creates a gamma-distributed delay in the activation of the proteins with mean time-delays B and C . In addition, noise in time-delays quantified by Coefficient of Variation where 〈. 〉 denotes the expected value. The activated protein > inhibits the expression of the protein > and the activated protein > enhances the production of protein G in the neighboring cell. It means that now the production rate of protein > is where L > ( ) and M G ( ) denote the count levels of L > and M G at time , respectively. Finally, the protein decays with rate . The overall model in Cell a consists of the following chemical reactions: (4) The model in Cell b is identical to that of Cell a.
We run Monte Carlo simulations to test our model numerically. Our code is written based on SSA algorithm (Gillespie, 1976) and we use the software StochKit2 to simulate the model (Sanft et al., 2011). In our code, we checked a wide range of time-delays. We start with timedelays B = 12.15 and C = 24.3 and we increased the time-delays up to B = 54 and C = 108 ., i.e. ≈ 4.5 fold change in time-delays as measured previously (Ay et al., 2014). We considered that = = 10, which means noise in time-delays is R S 4 = R W 4 = 0.1. Further note that, to maintain the oscillations of two cells in synchrony, we increase the time-delays with the same ratio. The protein degradation rate is = 0.2/ , which means the protein halflife is 3.46 as measured previously (Ay et al., 2013). Finally, we select the rates as B = 0.6 and C = 0.2, and by changing delay, we also change E to keep the mean of protein at 55 ± 1. Figure S5B illustrates that this model can replicate the synchronized oscillations by interplay between intracellular negative feedback and intercellular positive feedback. Our goal is to calculate the noise in a population of synchronized cells, however due to model size, increasing the number of cells results in increasing the complexity of model. Hence, instead of increasing the number of cells in our code, we used the peak values of oscillations in a cell as a proxy of a population of synchronized cells. To do so, after simulating the model using StochKit, we used a MATLAB code to process the data. In this code, we first read the data from the files that StochKit made. In the next step, we smoothened the data to remove the fluctuations contributed from random production and degradation events. The smoothened data was used to find the peak times. The difference between these times gave the period of oscillations.
Afterwards, we extracted the peak levels from the unfiltered data associated with these peak times. The expression variability at the peak levels was defined as expression noise. The