Synchronizing stochastic circadian oscillators in single cells of Neurospora crassa

The synchronization of stochastic coupled oscillators is a central problem in physics and an emerging problem in biology, particularly in the context of circadian rhythms. Most measurements on the biological clock are made at the macroscopic level of millions of cells. Here measurements are made on the oscillators in single cells of the model fungal system, Neurospora crassa, with droplet microfluidics and the use of a fluorescent recorder hooked up to a promoter on a clock controlled gene-2 (ccg-2). The oscillators of individual cells are stochastic with a period near 21 hours (h), and using a stochastic clock network ensemble fitted by Markov Chain Monte Carlo implemented on general-purpose graphical processing units (or GPGPUs) we estimated that >94% of the variation in ccg-2 expression was stochastic (as opposed to experimental error). To overcome this stochasticity at the macroscopic level, cells must synchronize their oscillators. Using a classic measure of similarity in cell trajectories within droplets, the intraclass correlation (ICC), the synchronization surface ICC is measured on >25,000 cells as a function of the number of neighboring cells within a droplet and of time. The synchronization surface provides evidence that cells communicate, and synchronization varies with genotype.

Supplementary Text S1 Experimental protocol for single cell measurements on the biological oscillator of Neurospora crassa • N.crassa strains, MFNC9 and MFNC30, growth 1.Inoculate a slant or multiple slants (depending on the growth characteristics of the desired strain) of growth media 1 with N. crassa conidia from a frozen culture.
2.Then put the cells at 30 °C in darkness to grow for about 48 h and then put them under light conditions on the bench top at room temperature. Keep the cells there for a minimum of 24 h, even 3-4 days to achieve maximum growth.
3.Then collect the conidia by suspension with sterile water (or other desired media). Determine the amount of suspension solution to be used by the growth and desired conidial concentration. 4.Filter the conidial suspension through a funnel containing a plug of sterile cotton into a sterile flask. Then pour the filtered conidial suspension into a 50 ml screw cap tube for storage. 5.Take a 20 µl sample of the conidial suspension and measure it on the cellometer (Nexcelom Cellometer Auto 2000, Nexcelom Bioscience LLC., Lawrence, MA) to determine the concentration of the conidia. (This measurement can also be taken with a haemocytometer under a light microscope.) 6.Then either dilute (with additional liquid) or concentrated (by centrifugation) the conidia to obtain the needed concentration.
• PDMS device preparation and surface treatment Fabricate the prototype polydimethylsiloxane (PDMS) microfluidic device for droplet generation through a standard soft-lithography approach and attach it to a glass slide 2 . The workflow of making the device is as follows.
2.Print the pattern of the photomask on the silicon wafer using a MA6 mask aligner (SUSS MicroTec Inc., Sunnyvale, CA) under a UV light source of wavelength 365 nm.
3.Develop the silicon wafer by submerging it in SU-8 developer (MicroChem Corp., Newton, MA) for several minutes. 4.Mix the pre-polymer Sylgard 184 (Dow Corning Corp., Midland, MI) with a cross-linker in a ratio of 10:1 w/w and then pour them onto the patterned silicon wafer and degas the mixture in a desiccator. 5.Put the mixture at 80 °C for 2 h for curing. 6.After curing, cut the PDMS layer and peel it off the wafer. Punch Inlets and outlets using a biopsy punch (Milte, Inc., York, PA) with an outer diameter of 1 mm. 7.Measure the thickness of the device by a profilometer (Dektak 150, Veeco Instruments Inc., Chadds Ford, PA). Normally, the thickness of the device is 50 µm.
8.Before attaching to a glass slide, treat the channel surface with air plasma (PDC-32G plasma cleaner, Harrick Plasma, Ithaca, NY) at 11.2 Pa O 2 partial pressure with 18 W power for 2 min. 9.After attachment, put the device in the oven at 80 °C for 8 h. 10.Afterwards put the device into a petri dish. Put the petri dish in a desiccator in a fume hood. Pipette 25 µL 1H,1H,2H,2H-Perfluorooctyltriethoxysilane (Sigma-Aldrich, St. Louis, MO) into the petri dish and let them surround the device.

•
Generate droplets with encapsulated cells and droplet collection 1.Dilute the original cell suspension to the needed concentration for encapsulating single-cells/multiple-cells in droplets. Use percoll as a dispersant to prevent the cells from settling down in the syringe during cell encapsulation 2.Insert ~2 cm of PE tubing into the two inlets of the device. Flush the device with DI water for the outlet to get rip of the PDMS debris.
3.Load the continuous phase into a 1 mL syringe. The continuous phase is a mixture of fluorinated oil FC-40 with 5.0 wt% of 008 fluorosurfactant. 4.Load the final cell suspension into a 1 mL syringe. 5.Connect the loaded syringes to the inlets of the devices with PE tubing and the tygon tubing. 6.Put the devices onto the stage of the microscope and secure the syringes to the slot of syringe pump. Set the flow rate of the continuous phase to be 13 uL/min and that of the dispersed phase to be 0.5 uL/min. and then run the syringe pump. 7.Connect the outlet of the device with a 1.5 mL conical tube using a PE tubing to collect the droplets in the conical tube. 8.After enough droplets are collected in the conical tube. Use a syringe to take out the redundant carrier oil (the fluorinated oil). Then use a capillary tube to retrieve the droplets for observation by submerging one end of the capillary tube in the droplets layer in the conical tube. 9.Put the capillary tube onto a clean glass slide and seal the tube with Epoxy glue. 10. To synchronize the cells, put the glass slide with capillary tube under the light source for 26 h.
• Image processing 1.Correction of CCD Imaging Imperfections 1.1 Dark current correction: Take ten images with the same exposure time as the timelapse experiment before each experiment but with the shutter closed. Average these images to create a dark current image. Subtract from each fluorescence image, the dark current image. 1.2 Bias correction: Take a few of images with zero exposure time and with the shutter closed. Average these images to create a bias image. Subtract from each fluorescence image (corrected for dark current), the bias image. 1.3Flat field correction: Take a couple of images with quite short exposure time of uniformly illuminated screen/field of view. A Rhodamine B solution with a concentration of 6.6 µM in a 35 mm culture dish was used as the uniformly illuminated screen. Average these images to create a flat field image. Obtain the average pixel value within the (dark/bias corrected) flat field image (call it α). Then, for each fluorescence image (dark/bias corrected), multiply all pixels by α and then divide (pixelwise) by the flat field image. Flat field image should be created each time optical properties of the illumination and lens system are changed. For example, a different flat field image is used for different magnifications.

2.Fluorescence intensity data extraction
We wrote a Matlab routine to sort droplets according to the number of cells in them and to track the fluorescence intensity of individual cell frame by frame. Locations of cells are determined by finding a local maximum in fluorescence intensity. The sorting is based on the difference of distance between cells and droplet and the radius of the droplet. The tracking is mainly based on the minimum distance between cells in consecutive frames. The quality control filter to filter out cells that are not linked property is based on comparing values of the extracted fluorescence intensity in consecutive time for one cell and between cells. The quality control filter marks the cells that are not linked properly as bad cells and they are not considered later in the analysis.
Additional controls recommended are an examination of a time lapse video S1 for detecting germination, fusion of cells, and cell division. In addition we selected ~200 singletons and monitored surface area over time for germination and cell division (10 are shown in Fig. S5). In addition, we plotted ~200 random pairs of cells and tracked them over time to see if there was any change in distance between them that might be evidence of cell fusion (10 are shown in Fig. S5).
We also recommend an examination of a phase histogram of single cells (Fig. 3) to validate that the phase histogram is unimodal. Bimodality might arise if cells were not all synchronized at the same time. A dip test was performed (P = 2.2 x 10 -16 ) 3 . The fact that a tiny dip of the same size also occurs in the stochastic simulation (Fig. 3D) may suggest that the dip test does not capture the intracellular stochastic variation in the data and that we may need to speed up encapsulation to less than 20 m, the current time taken for the encapsulation process. A third control is to substitute doped beads for cells and observe the beads over 10 days in the capillary tube. The periodogram of single beads should be computed (Fig. 3). We found a maximum in the periodogram at 20 h of about 0.03 (more than 4-fold less than the maximum in periodogram in Fig 2D), which is explainable by the 1 o C fluctuation in temperature occurring diurnally during the microfluidics experiment. A fourth control is to examine the period of ensembles of oscillators in various strains in race tubes to confirm that the period of each ensemble is the same as that of individual cells measured by microfluidics. A final control is to examine race tubes exposed to a flash on a iPhone every 30 minutes at the same intensity as the exposure of single cells to light every 30 minutes to observe fluorescence. We confirmed no change in the clock phase and period in race tubes relative to tubes without a flash exposure every 30 m. The banding with the light show some irregularity compared with those race tubes in total darkness, but the light exposure was longer (~0.5 seconds) in Fig. S6.
S2 MCMC procedure developed for fitting a stochastic genetic network

•
Stochastic Model of the Clock's Intra-Cellular Reaction Kinetics The time evolution of a well-stirred chemical system, such as the biological clock, is usually described by a system of ordinary differential equations (ODEs) that show the evolution of molecular numbers of species in the system as a continuous, deterministic process 4 . The problem with this approach is that molecular numbers are given by whole numbers, and when they change, they always do so by discrete amounts. If molecular numbers of some species are small or if the dynamics of the chemical system is susceptible to noise 5 , then a system of ODEs is unlikely to describe completely the dynamics of such a system . We need a discrete, stochastic model to describe the noisy behavior of the system, and we need a stochastic model to describe how chemical reactions occur as discrete events at the molecular level. Such a model, and an algorithm to solve it, was devised by Gillespie 6 .

•
Gillespie algorithm We consider a system of chemical species interacting through chemical reactions . The system is assumed to be confined to a constant volume. The number of molecules of species at time , denoted by , is an integer-valued, piecewise constant function of time. undergoes discontinuous jumps at random points in time, whenever a reaction event occurs which consumes or produces molecule(s) of species . To simulate the evolution of the state vector, , we need to be able to calculate two things: how much time will pass until the next reaction event occurs and which reaction will occur at that time.
Each reaction is characterized by two quantities. The first is its state-change vector , where is the change in the molecular number of species , caused by one -reaction event. These state-change vectors are defined by the topology of the reaction network 4 , as summarized graphically in Fig. 4A. If the system is in state at time and reaction occurs, then the system instantaneously changes to state . The second quantity characterizing is its propensity function Given a state , is the probability that one reaction occurs in the next infinitesimal time interval The following assumptions, based on physical principles, allow us to determine the distribution of the next reaction time and the probability that is the next reaction.
(1) =probability that reaction takes place in the infinitesimal time interval .
(2) , where =reaction rate coefficient of and =number of possible ways can occur given system state . (This is a mass action assumption used in the working model of the clock now supported by over 30,000 data points 7 ). The rate constants k j are listed in Fig. 4A under the no coupling assumption.
For instance, is , and for zeroth-order, unimolecular, hetero-bimolecular and homo-bimolecular reactions, respectively, involving species and/or for given the state prior to the reaction event.
For a given model parameter set, i.e., given rate coefficients and initial conditions, we generate over 1000 random trajectories of the stochastic clock model using the Gillespie Algorithm 6 . Each trajectory starts from the given set of initial molecular numbers for all species, , at time . Random-length time steps, from one reaction event to the next, are then performed, based upon foregoing assumption (1), that the time until the next reaction, , is exponentially distributed with a rate (S1) where t is either the time of the most recent reaction event, or else . Thus, a random is drawn from that exponential distribution. The probability that this next reaction is , is given by Thus, to decide which reaction, , takes place next, we choose a random number between 1 and N R according to these probabilities. Namely, we draw a uniform random number between and and find that index value for which . (S3) Then, the reaction that occurs next is . Time is then advanced to that next reaction event by a time step of length .The trajectory is terminated when the last reaction event time, , exceeds the total model simulation time duration, . In summary, the Gillespie trajectory generation algorithm works as follows: typically with . Here, denotes the rate coefficient of reaction and is the (generally unknown) unit conversion factor from the fluorescent protein molecule count into the corresponding normalized and de-trended fluorescent photon count signal, as measured experimentally by the CCD and explained below. The mean is then taken over these simulated trajectories, e.g., for the mean bare periodogram, as follows: and the periodogram (power spectrum), of the th Gillespie trajectory is given by Here, X n obs ,k (t j ;Θ) is the molecule count of the observable (i. • Fluorescent Signal Data Normalization Several corrective data transformations and a detailed noise analysis must be performed with the observed single-cell fluorescence time series raw data, before they can be Fourier analyzed to obtain periodograms which can then be meaningfully compared to the stochastic model predictions. The basic detection method for time-varying single-cell protein concentrations is to record the fluorescent single-cell signal, generated by the protein molecules upon illuminating the cell with radiation from a shorter-wavelength light source. Let denote the mean fluorescent photon count signal produced by The proportionality constant is a property of the fluorescent protein molecule, i.e., depends only on the fluorescence properties of a single fluorescent protein molecule in its intra-cellular environment, but is assumed to be the same across all cells at all observation times. The time and cell indices, respectively, are j = 1,..., L and . Here, is again the total number of equidistant observation time points t j = ( j −1)T / (L −1) , and is the total number of cells. For the experimental analysis and model ensemble simulations reported here we used T =85h, , and , unless stated otherwise, corresponding to an observation time interval t j+1 − t j = 0.5h.
Since the illumination intensity produced by the lamp can vary with time over the 200+ hour duration of the experiment, a reference measurement is made, almost at the same time as each fluorescent cell observation, to also detect the fluorescent intensity, , from a Rhodamine B (RB) sample illuminated with the same lamp, and hence also proportional to the illumination photon count, : , is a property solely of the RB sample and assumed to be constant in time. Hence, the time-dependence of the single-cell signal, , due to the time-varying illumination can be cancelled out by normalizing the single-cell signal, , to the reference intensity, , i.e., by taking (S8) which is proportional only to the fluorescent protein count in cell # at time . The foregoing results now need to be extended to allow for the presence of detection noise.
Note in passing that the ratio a C / a R in Eq.(S8) is related to the periodogram unit conversion factor φ in Eq. (S5) by: φ = (a C / a R ) 2 . Since that ratio is unknown, so is φ . Hence, φ is treated as an unknown, adjustable model parameter, along with all other stochastic model parameter variables in . •

Intra-Cell Stochasticity and Detection Noise in Observed Signals
The cellular fluorescent signals, , in Eq. (S6) are random variables due to the stochasticity of the intra-cellular chemical processes. They reflect the stochastic nature of the underlying chemical reaction netwerk, via the fluorescent protein molecule counts, , in Eq. (S6). By definition, , does not contain any randomness from the experimental detection process. The detection process adds further randomness in producing the CCD detector output, denoted by and consisting of the stochastic detector input signal, , superimposed with a detection noise term, : . (S9) In the following, all cell-related fluorescence signal and respective periodogram quantities containing random contributions from experimental detection noise will be Corresponding quantities that would be observed (hypothetically) in the absence of any experimental detection noise are denoted without the tilde overscript, such as The latter quantities exhibit randomness only due to the stochasticity of the intra-cellular chemical processes, as reflected in .
Each detection noise variable, , in Eq. (S9) is assumed to be statistically independent across all cells, , and all times , with zero mean and a detection noise variance : (S10) where and are Kronecker delta-functions. The bracket denotes the conditional population average over all possible experimental detection noise configurations, given the "clean" detector input signal, , prior to the addition of detection noise. That is, is an average over a (hypothetical) infinite number of repeats of the detection process, each repeat detection performed with the same, fixed detector input, . Thus, Eq. (S10) is not averaged over, but conditional upon, the stochastic intra-cellular state trajectories, , and the resulting noise variances, , depend on the stochastic detector input signal , given by Eq. (S6).
The detection noise likewise modifies the illumination-normalized pre-detection signal, , with a noise term denoted by , resulting in an illumination-normalized postdetection signal , (S11) again with statistically independent noise variables, ε k t j We are neglecting the detection noise effect on the reference signal, . The latter is recorded by a large number of pixels, across the entire field of view of the detecting CCD camera and it is also not subject to variability from random focal point excursions. Based on the detection noise analysis described above and below, the reference signal can thus be assumed to have a much lower detection noise level, as a percentage of signal, than the cell-generated detector output signals .
Intra-cellular stochasticity and the detection noise both contribute, statistically independently, to the random variability of the actually observed CCD detector output signal . It is necessary to separate and account for both sources of randomness in our modeling, since only the intra-celluar stochasticity, but not the detection noise, is of biological interest. As discussed below, the detection noise also introduces a bias into the periodograms constructed from the observed output signals, . A separate estimate for the detection noise variance is thus needed to correct for this bias.
There are six sources of detection randomness in single cell measurements that are associated with the fluorescent signal generation and detection from a cell, as follows: (1) the quantum randomness of fluorescent photon generation in the cell; (2) the cell's random focal plane excursions; (3) the quantum randomness of photo-electron generation in the CCD camera; (4) CCD random pixel-to-pixel variations in each pixel's photon-toelectron conversion yield; (5) CCD image pixel summation; (6) CCD bias and dark current subtraction. For the analysis that follows, these six noise sources are assumed to be statistically independent and to add noise successively, with the noisy output signal from one noise source serving as the input signal of the next. Furthermore, we assume, on simple physical grounds, that each of the foregoing six noise sources can be classified as either 0 th , 1 st or 2 nd order, in the sense that the variance of the noise added to the input signal by a 0 th , 1 st or 2 nd order source is proportional to the 0 th , 1 st or 2 nd power of the source's input signal, respectively. Of the foregoing six noise sources, source (6) is 0 th order, sources (1) and (3) are 1 st order, and sources (2), (4) and (5) are 2 nd order. A detailed mathematical noise compounding and propagation analysis then shows that the successive application of any sequence of 0 th , 1 st or 2 nd order detection noise sources to the detector input signals can all be lumped into a single "effective" noise source, with statistically independent "effective" noise variables, , whose respective variances depend on the original detection input signal by a simple quadratic relation of the form: . (S13) The coefficients of this quadratic function, , and , can be expressed in terms of the various proportionality coefficients of the underlying individual 0 th , 1 st or 2 nd order sources. These proportionality coefficients can not be estimated from available experiments for each individual source. , and must therefore be determined by a separate "noise calibration" experiment where the intrinsic noise variance of the detector input signal is either zero, or can be quantified separately by other means. The analysis and results of these noise calibration experiments, using samples of fluorescent beads instead of living cells, will be described in the following section below.
To simplify the periodogram noise analysis discussed below, we make the additional assumption that the variances can be approximated by their population mean over the stochastic intra-cellular molecule count time series, X k (t j ) , as measured by the mean photon count signal in Eq. (S6), i.e., by Eq. (S13): where ... c denotes the population mean over the intra-cellular stochastic time series, The resulting means S k (t j ) c and (S k (t j )) 2 c , and hence (σ ξ, j ) 2 , are independent of the cell index k, since, by model assumption, the random distribution of the stochastic intra-cellular time series, X k (t j ) , is the same for all cells.
From Eqs. (S9), (S10) and (S13), we then get denotes the joint population average over detection noise configuration and over intra-cellular stochastic time series variable, X k (t j ) . For large cell sample sizes, K, the foregoing population averages can be approximated by sample averages over the actually observed CCD signals from all cells: Eqs. (S16) and (S15) were then used to estimate the detection noise variances (σ ξ,k, j ) 2 .
• Detection Noise Calibration by Bead Experiments To quantify the detection noise level in the fluorescent single-cell signals by way of the quadratic signal-to-noise relation, Eq. (S15), we performed a series of noise calibration experiments where the sample of droplet-encapsulated living cells was replaced by a nearly mono-disperse sample of small, droplet-encapsulated, fluorescent beads. The beads were polymer microspheres internally doped with fluorescent dye, Cat. Code: FS06F, Envy Green (Excitation/Emission:525nm/565nm), with mean diameter 9.94 μm and diameter standard deviation 0.76 μm, as provided by the manufacturer. In contrast to the cellular fluorescent signals, the mean fluorescent output signal from each bead is not subject to any intrinsic time-dependent stochastic variability. However, the beads do exhibit time-independent variability in their fluorescent signal output due to a small, known, but non-negligible variance in the bead size distribution. The overall variability in the CCD detector output signal observed in the bead experiment is thus the result of the experimental detection noise and of the bead size variability. The latter can be treated as an added 2 nd order noise source, while the six noise sources from the detection process are the same as those in the single-cell experiment. By subtracting out the signal variability caused by the random bead sizes, we can thus extract from the bead experiments the noise calibration coefficients, , and , entering into Eqs. (S13)-(S15), as follows: Let and denote, respectively, the mean detector input signal and the recorded CCD detector output signal for the th bead in the sample, both at observation time . (S17) The mean bead signal, , represents the mean fluorescent photon count per bead emitted into the CCD detector at observation time , averaged over all possible bead sizes. It therefore does not depend on the bead index . The bead noise variables, , contain the random contributions from both the detection noise and from the bead size variability. They are assumed to be statistically independent across all beads (but not across all times!), with zero mean and a total noise variance : Here, the brackets denote the joint population average over all possible experimental detection noise configurations and over all possible bead sizes.
Applying the same noise compounding and propagation analysis as for Eq. (S13), but including both the same six detection noise sources as in Eq. (S13) and in addition the bead size variablility, we can then show that (S19) where coefficients , and , in Eq. (S13) are related to , and , , by , , .
Here, is the relative variance of the emitted fluorescent photon count per bead, due to bead size variability only, in the absence of any detection noise. This is discussed in more detail below. From the mathematical noise compounding and propagation analysis it follows that the coefficients , and , in Eq. (S13), as well as , and , in Eq. (S19) must be non-negative. This requirement will allow us to impose important constraints on the bead parameters, as discussed below.
For a sufficiently large bead sample size, , we can then estimate the -values and corresponding -values from sample means of the observed bead data, : where is the total number of beads observed in our experiment at time . The fluorescent photon counts in the bead experiments bracketed those used in the cell experiments. Also shown in Fig. S7 is the result of a least-squares fit of Eq. (S19), subject to the non-negativity constraint imposed on each fit parameter, , and . Fig.  S8 shows an analogous plot and fit for the data from Experiments 1-5 only. Both fits were remarkably close to the data, both with R 2 > 0.99. The fit parameter results are summarized in Table S1. It is important to notice that the -term dominates in Eq.
(S19), relative to the -and -terms, over the relevant range of -values, which is of order 10 6 , measured in our CCD camera output units.
To estimate the parameter , we need to model how each bead's emitted photon count, , depends on the bead size, quantified in terms of the bead diameter, . This dependence is controlled by the spatial distribution of the fluorescent dopant inside the bead. As a simplest model of the bead doping profile, we assume that some fraction of the dopant material is uniformly distributed, with a -independent bulk dopant density, over the 3D spherical volume of the bead, while the remaining fraction of dopant material is uniformly distributed over the bead's 2D spherical surface, with -independent surface dopant density. The total number of fluorescent photons, , emitted by a bead into the CCD detector, is assumed to be proportional to the total dopant amount in the bead and proportional to the number of incident illumination photons, , during exposure. Hence, we can write as a function of the form . (S23) The first and second term in Eq. (S20) represent, respectively, the volume and surface dopant contributions to the total fluorescence: they are proportional to the bead volume, , and to the bead surface area, , respectively, rescaled by the  Table S1 for the two -values extracted from the fits of Eq. (S19) to the Experiment 1-5 and Experiment 0-5 bead data sets. As discussed below, other non-negativity constraints will also impose very tight lower bounds, and σ B,min 2 , which then limits and σ B 2 to an even narrower range of allowed values, as also shown in Table S1. We therefore adopt the "volume fraction" variable, , as a convenient model input parameter to quantify the effect of the bead doping profile on the variability of the bead fluorescent photon output, via Eqs. (S26) and (S20).
• De-Trending of RB-Normalized Time Series Data A first step in any periodogram analysis is to ensure that the underlying time series of fluorescence for each cell is stationary. Because of photobleaching (Fig. 2B), the time series of ccg-2 fluorescence for each cell was not stationary. It was therefore necessary to carry out a moving average de-trending 8 on the RB-normalized cell time series data, , from Eq. (S11). The first step in this procedure was to subtract from the data a linear trend line, , fitted to the data over short contiguous time blocks, each block containing successive time points, , as follows: otherwise. Here, the block index, , runs over the values . The regression coefficients, and , in this moving average de-trending were estimated by linear least-squares fit to the data points from the respective time block, : and The final de-trended signal, denoted by , was then obtained by averaging the block- , from all those time blocks which contain the time point , i.e., over all -values for which . Combining that averaging procedure with Eqs.
(S28), (S29) and (S30), the resulting can be written as a linear transformation of the non-de-trended RB-normalized CCD output data, , with a transformation "weight" matrix as follows: and The th block indicator function, , is given by and the mean block time is The quantity is the number of time blocks, , which contain observation time : By Eq. (S11) we then get, for cell # k at observation time t j : Here, ζ k (t j ) is the de-trended, RB-normalized experimental detection noise and V k (t j ) is the de-trended, RB-normalized fluorescent signal in the absence of experimental detection noise: The noise variables ζ k (t j ) are not statistically independent anymore at different times and their (co-) variance matrix elements, obtained via Eqs. (S43) and (S12), are The variances (σ ξ, j ) 2 in Eq. (S48) were estimated from the observed CCD data by Eq.

Observed Cell-Averaged Periodogram and Periodogram Variances
The cell-averaged periodogram in Eq. (4), Q obs ( f  ) , that can be directly compared to the model prediction, Q exp ( f  ) , must be defined in terms of the de-trended, RB-normalized signals, V k (t j ) , that would be observed in the absence of detection noise, as follows: Here, ... c again denotes the population mean over the intra-cellular stochastic time series, X k (t j ) and and is the periodogram of the fluorescent signal V k (t j ) and F k ( f  ) is its Fourier transform, at frequencies f  =  / T with  = 0,1,..., L and imaginary unit i with i 2 ≡ −1 : Note that, in Eq. (S49), the population mean, Q k ( f  ) c , taken over all intra-cellular stochastic time series, X k (t j ) , is again independent of the cell index k, analogous to Eq. (S14).
The de-trended signal,V k (t j ) , in Eq. (S52) is based, by linear transformation, on the "clean" fluorescent signal, Y k (t j ) , without experimental detection noise. Therefore, Q obs ( f  ) , cannot be directly calculated from its defining Eqs. (S49)-(S52). Instead, we must express Q obs ( f  ) in terms of the actually observed, RB-normalized CCD signal, , which is subject to experimental detection noise. The result is as follows: Here, ... e,c again denotes the joint population mean over both the experimental detection noise and over the intra-cellular stochastic time series, X k (t j ) , and with  V k (t j ) calculated by Eq. (S11) and Eqs. (S35)-(S40), from the observed CCD in Eq. (S54) does not depend on the cell index k.
(S65) is negligible compared to theσ ε 2 -term. From our experimental data, we find that σ ε,max1 2 is the most restrictive, i.e., the lowest upper bound onσ ε 2 , and hence, it is the overall upper bound:    Table  S1. Notice that Q obs ( f  ) -values become negative as f B is lowered below f B,min .
However, initially, for f B -values below, but still close to f B,min , this happens only at the shortest periods ( ). As f B is lowered further, the range of negative-Q obs ( f  ) periods with increases.
The top two Q obs ( f  ) -curves in Fig. S9 show what happens when f B is varied over the allowed range: from f B,min ≅ 0.7679 to f B,max ≅ 0.7775, Q obs ( f  ) varies by no more than 10% of its peak value, max  Q obs ( f  ) ! " # $ ≈ 0.04, across all periods T  ≡ 1 / f  . Clearly, Q obs ( f  ) is not substantially affected by the choice of f B -value within the allowed range T  ≡ 1 / f  As an alternative approach to estimating σ B 2 and f B we also exploited the time dependence of the relative bead noise correlations, defined by Thus, s (b) ( j, j ') is the temporal correlation between the noise variables η k (b) (t j (b) ) and , on the same bead, k = k ' , at bead observation times t j [Note here that By averaging over all beads in the six bead data sets, Experiments 0-5 and replacing population averages ... e,b in Eq. (S73) by sample averages analogous to Eqs. (S21) and (S22), we are able to extract reasonable estimates for .
The basic idea of the σ B 2 -estimation from time correlations then rests on the observation that, on physical grounds, the correlation contribution to s (b) ( j, j ') from all the six detection noise sources will die out as a function of temporal distance, since the detection noise sources are subject to random temporal fluctuations. On the other hand, the "noise" contribution from the bead size variation does not die out, since the bead size, D , does not fluctuate with time. As a simple physical model, we assume that the detection noise correlations die out exponentially with Δt( j, j ') , i.e., with a functional dependence of the form The constants a s , τ s and b s denote, respectively, the equal- contribution from the temporally fluctuating detection noise sources, the exponential decay lifetime of the temporally fluctuating detection noise, and the contribution from the temporally constant bead size variation. In the limit Δt( j, j ') → ∞ , only the bead size variance contribution to survives. Recall now that, is the relative variance of the emitted fluorescent photon count per bead, due to bead size variability only, in the absence of detection noise. The model constant b s is then directly related to the relative bead variance parameter, σ B 2 , by We therefore extracted s (b) ( j, j ') from the six bead data sets obtained in Experiments 0-5.
We then fitted those s (b) ( j, j ') -data with a fitting function of the form given by Eq. (S75), with fit parameters, a s , τ s and b s . The resulting values for σ B 2 from Eq. (S76), were σ B 2 = 0.04778 ± 9.3×10 −5 if all six Experiment 0-5 data were included, and σ B 2 = 0.03864 ± 1.1×10 −4 if only the five Experiment 1-5 data were included,. The given standard deviations of were estimated by bootstrapping the respective bead data sets over all beads, k = 1,..., K (b) , for each of the six or five bead experiments, and then performing a separate fit on every bootstrap data set thus generated. As seen in Table S1, the resulting two estimates bracketed the desired range σ B,min and they were within 7% of that range. We conclude that the basic model assumptions underlying our detection and bead size noise analysis is well supported.
was used, with Ω denoting the normalization factor. The objective of the ensemble simulation is then to generate a large random sample of parameter variable vectors, Θ , drawn from , which provide near-optimal fits to the experimental data.
In principle, we could have constructed from the squared residuals of the biascorrected bare (un-normalized) observed periodogram, from Eq. (S53), and the corresponding bare expected periodogram, Q exp ( f  ;Θ) from Eq. (S4). However, in light of the relative smallness of the above-estimated detection noise contributions to the variance and bias of the observed bare periodogram, we employed a simplified approach which neglects these small detection noise corrections. This simplified approach also eliminates the uncertainties due to the unknown unit conversion factor variable φ entering into the model predictions via Eq. (S5). This is achieved by defining in terms normalized instead of bare periodograms, as was done in Eq. (4). The normalized observed periodogram in Eq. (4) is defined by Likewise, the normalized expected periodogram in Eq. (4) is defined by where Q k exp ( f  ;Θ) is given by Eq. (S5) and the normalization factors are Lastly, the experimental variances in Eq. (4) are given by (S82) Since i.e., The Q obs ( f  ) and Q exp ( f  ;Θ) can therefore be compared directly, without any unit conversions, as was done in Eq.(4).
Standard Metropolis MC updating steps, based on L(Θ) , were performed on randomly selected single components of the parameter variable vector, , At the start of the simulation, the X n (0) -and k m -variables in Θ were initialized by appropriately rescaled initial concentration and rate coefficient values, respectively, obtained from earlier ensemble simulations of the corresponding deterministic model. 10 (  Fig. 2D is the fit corresponding to the minimal = 130.231. are the best fit of Eq. (S19), 2 , subject to the three positivity constraints The three positivity constraints effectively introduce non-linearity onto the least-squares fitting procedure and can result in best fit parameter choices where one or two of the three fit parameters are forced to assume values of zero exactly, as found in all three fits reported here. The fitted input data, were estimated by Eqs. (S21) and (S22) from the bead fluorescence intensities observed in Expts 0-5. Standard deviations (Std.Dev.) were computed by bootstrapping the bead experiments and repeating the fit for each bootstrap sample. From B (b) fit , we estimated (σ B,,max ) 2 by Eq. (S30). From Eq. (S72), we estimated (σε ,,max ) 2 , the upper limit on σε 2 , which was then used to impose the upper limit on the B-coefficient in Eq.(S15), denoted by B max , via Eqs. (S15), (S48) and (S47). Eq. (S20) and B max were then used to impose the lower σ B 2 -limit, denoted by (σ B,min ) 2 . The values of f B,max and f B,min were then calculated from (σ B,,max ) 2 and (σ B,min ) 2 , respectively, by Eq. (S26). The best fit [1], to the six-experiment bead data set, Expts 0-5, resulted in an empty allowed σ B 2 -range , i.e., (σ B,min ) 2 > (σ B,max ) 2 or, equivalently, f B,min >f B,max . The (σ B,min ) 2 and f B,min from fit [1] also exceeded their physical upper limits, (σ 2 B,Hi ) 2 0.052613 and 1, respectively, imposed by Eq. (S27) and (S24). An additional, more restricted fit [2] to the six-experiment bead data set, Expts 0-5, was then performed, with A (b) forced to A (b) =0, and resulted in a finite allowed σ B 2 -range , i.e., (σ B,min ) 2 < (σ B,max ) 2 . The best fit [3], to the five-experiment bead data set, Expt 0-5, resulted in a finite allowed σ B 2 -interval , i.e., (σ B,min ) 2 < (σ B,max ) 2 , without additional restrictions on the fit parameters beyond the three positivity constraints stated above. The results from best fit [3] were used for the subsequent detection error and bias analysis of the single-cell periodogram data.  S1.
Histograms of the times of activation are measured in a window from 30 hours to 115 hours for 1,024 Gillespie trajectories under the best fitting stochastic model, χ 2 = 130.231. The times of activation of these genes are the times of the reaction frq 0 à frq 1 or ccg-2 0 -> ccg-2 1 -30 hours in each of K=1,024 realizations of the Gillespie stochastic trajectories generated by the model with minimal = 130.231 using the Gillespie Algorithm described in Eqns (S1)-(S3). The histogram counts of these two reactions are of all occurrences of frq 0 à frq 1 (Panel A) and of ccg-2 0 -> ccg-2 1 (Panel B). All times of activation are measured from the 30 hour starting time in all 1,024 Gillespie Trajectories. The times of the reactions are only recorded on the time interval from 30 hours to 115 hours so that the time of gene activation can vary from 0 to 85 hours. The zero point in the graph corresponds to 30 hours, and maximum time of activation recorded corresponds to 115 hours.

Fig. S2
The synchronization as measured by the "order parameter" R in Garcia-Ojalvo 10 or Hilbert transform phase synchronization 11 S3. The Kuramoto order parameter K displays an alternating structure with number of cells per droplet and substantial coherence between oscillators. The order parameter K is the mean absolute deviation of oscillators phases 12 . The last 30 hours of a ten day microfluidics experiment were discarded, and all droplets with less that 80% successful calls on tracking over time were not used. The phase was computed about zero by removing the mean, extracting the discrete Hilbert Phase, and continuizing the discrete Hilbert phase as described in the legend of Fig. 3.    Table S1. The best fit parameter values, A (b) , B (b) and C (b) , are recorded in Table S1 as Fit [1].
.  Table S1 as Fit [3].  Table S1, Fit [3]. Each curve is a periodogram, Q obs ( f  ) , corrected for bias due to detector noise according to Eqs. (S53)-(S61), assuming different values of bead volume fluorescence fraction, f B and the noise calibration coefficients A (b) , B (b) and C (b) recorded in Table S1 as Fit [3].