article: Stochastic noise at the cellular level has been shown to play a fundamental role in circadian oscillations, influencing how groups of cells entrain to external cues and likely serving as the mechanism by which cell-autonomous rhythms are generated. Despite this importance, few studies have investigated how clock perturbations affect stochastic noise—even as increasing numbers of high-throughput screens categorize how gene knockdowns or small molecules can change clock period and amplitude. This absence is likely due to the difficulty associated with measuring cell-autonomous stochastic noise directly, which currently requires the careful collection and processing of single-cell data. In this study, we show that the damping rate of population-level bioluminescence recordings can serve as an accurate measure of overall stochastic noise, and one that can be applied to future and existing high-throughput circadian screens. Using cell-autonomous fibroblast data, we first show directly that higher noise at the single-cell results in faster damping at the population level. Next, we show that the damping rate of cultured cells can be changed in a dose-dependent fashion by small molecule modulators, and confirm that such a change can be explained by single-cell noise using a mathematical model. We further demonstrate the insights that can be gained by applying our method to a genome-wide siRNA screen, revealing that stochastic noise is altered independently from period, amplitude, and phase. Finally, we hypothesize that the unperturbed clock is highly optimized for robust rhythms, as very few gene perturbations are capable of simultaneously increasing amplitude and lowering stochastic noise. Ultimately, this study demonstrates the importance of considering the effect of circadian perturbations on stochastic noise, particularly with regard to the development of small-molecule circadian therapeutics.
  Circadian rhythms are daily changes in gene expression and physiology that persist even in the absence of external environmental cues [1]. In mammals, such rhythms are organized in a hierarchical fashion: at the tissue-level, the brain’s suprachiasmatic nucleus (SCN) serves as the master pacemaker and keeps circadian oscillations in peripheral tissues in phase with the light-dark cycle. In the SCN, cell-to-cell coupling keeps individual cells in tight synchrony [2], while coupling between circadian oscillations in peripheral tissues in vivo or cultured reporter cells in vitro is thought to be very weak or absent entirely [3,4]. Within each tissue, cellular-level rhythms in gene transcription are generated by a large network of interacting gene regulatory elements, in which time-delayed transcription-translation negative feedback gives rise to sustained oscillations [5]. The robust oscillation of circadian factors has been linked to metabolic health [6], since rhythms compromised by gene knockout [7] or irregular feeding schedules [8] result in an increased risk of metabolic disease. Additionally, as the amplitude of circadian transcription can be affected by lifestyle variables such as diet, age, or work schedule, there has been recent interest in developing pharmacological strategies for increasing the amplitude of circadian cycles in metabolic tissues [9]. A detailed understanding of the underlying transcriptional mechanisms is essential for the development of circadian therapeutics to be successful. The functional roles of different genes in circadian regulation have traditionally been studied using behavioral-level data and genetic knockout experiments [10]. Bioluminescence-based cellular circadian reporters offer a more direct view of the gene regulatory network [11] and are amenable to high-throughput screens, allowing genome-wide exploration into factors that affect circadian rhythmicity [12]. Additionally, cultured circadian reporter cells allow the change in transcriptional amplitude following a perturbation to be quantified. This additional parameter has proven useful in differentiating between perturbations with the same effect on period [13] and has aided the search for small-molecule therapeutics to boost clock amplitude [9]. Bioluminescence rhythms at the cell culture or tissue-level are the result of the collective behavior of thousands of individual cells. Transcription at the single-cell level is strongly affected by intrinsic cellular noise, caused by the low molecular counts of the mRNA and protein species involved. As a result, bioluminescence traces of individual cells are stochastic, with significant variability in both amplitude and period length from cycle to cycle [14]. In addition to intrinsic noise, circadian oscillations are also affected by extrinsic noise sources. Extrinsic noise results from heterogeneity between cells, such as differences in size or physical environment, leading to differences on a cell-to-cell basis in expected period and amplitude. The effects of noise in biological systems has been well-studied, and relative amounts of intrinsic and extrinsic noise can be identified through carefully designed single-cell experiments [15]. For circadian systems, intrinsic noise has been suspected to play a larger role: a single cell’s variability in period from cycle-to-cycle is larger than the variability in mean period length between cells [2]. However, both sources of noise have an effect on population-level rhythms: in cell cultures that lack cell-to-cell coupling, it has been shown that stochastic noise is manifested in damped oscillations at the population-level as individual oscillators gradually drift out of phase [14,16]. This type of behavior has also been seen in other experimental systems, such as NF-κB signaling or yeast glycolytic oscillations [17,18]. The amount of noise in system is therefore linked to the ability of tissue-level clocks to maintain high-amplitude rhythms. Despite the averaging that occurs at the population-level, cell-autonomous stochastic noise plays an important role in determining the function of the circadian oscillator. Noise in circadian rhythms has long been considered an important factor in how circadian rhythms have evolved [19]. A recent study showed that stochasticity is critical to the population-level response to a neuropeptide and forms the basis for how the SCN entrains to light-mediated cues [20]. Additional studies have suggested that the basis of single-cell rhythmicity may depend on stochastic noise, as models of deterministically damped oscillators, when simulated stochastically, capture the noise characteristics seen in single-cell fibroblast data equally well as limit-cycle oscillators [21]. Despite the importance of single-cell stochasticity in circadian rhythms, measuring stochastic noise currently requires careful preparation, recording, and image processing of individual cells [22]. As a result, while circadian perturbations have been postulated to affect single-cell stochasticity [23], no study has experimentally quantified changes to stochastic noise as a result of a small molecule or genetic perturbation. In this study, we demonstrate that changes to stochastic noise can be reliably inferred from the changes to the damping rate of population-level bioluminescence recordings of cultured circadian reporters. Our method assumes that oscillations in individual cells are both which have been shown to hold experimentally for cultured fibroblast cells. We demonstrate the validity and usefulness of such an approach on several types of circadian data. First, we show using single-cell fibroblast data that intrinsic stochastic noise is directly related to population-level damping. Next, we show that a small-molecule modulator is able to change damping rate in a dose-dependent fashion, and verify using a mathematical model that changes to intrinsic stochastic noise is a likely mechanism. Finally, we calculate the genome-wide effects of siRNA knockdown on overall stochastic noise, and demonstrate that population-level damping rate is independent of other circadian parameters, such as period or amplitude. Using this additional information, we show that circadian rhythms have likely evolved to an optimal point of high amplitude and low stochastic noise. Our results should prove especially important in the future search for small molecule circadian therapeutics, as it allows the effect of candidate drugs on stochastic noise to be quantified in a high-throughput manner.
  Raw experimental data x (ti), i ∈ {0, …, N − 1} are first detrended using Hodrick-Prescott filter with a smoothing parameter γ = 0. 05 (24hrs T s) 4, in which Ts is the sampling rate (in hours) [24]. The detrended data are then filtered using a low-pass filter to remove high-frequency noise (forward-backward Butterworth filter with n = 5, wc = 0. 1). We denote the detrended and filtered experimental data by y (ti). A damped sinusoid, specified by: y ^ (t) = A e - d t sin (2 π t T + θ) is then fit to the filtered data. For numerical efficiency, the period, T, and damping rate, d, parameters are fit first using a matrix pencil method [25], reviewed in [26]. Amplitude, A, and phase, θ, parameters are subsequently fit using a linear least-squares regression. Note that in this manuscript we use amplitude to denote the initial rhythm strength, and damping rate to denote the rate at which this strength decays with time. Overall R2 values for the regression were calculated from the residual error between the detrended data and fitted sinusoid: R 2 = 1 - ∑ i = 0 N - 1 (y (t i) - y ^ (t i) ) 2 ∑ i = 0 N - 1 (y (t i) - y ¯ (t i) ) 2 in which y ‾ (t i) represents the mean of the detrended data. Single-cell bioluminescence data for 79 cells was obtained from Leise et al. , 2012 [22]. As was done in the original study, a discrete wavelet transform (using PyWavelets, http: //www. pybytes. com/pywavelets) was performed to detrend and remove noise. A discrete wavelet transform decomposes the signal into multiple frequency bands [27]. By only considering frequency bands close to the circadian frequency, high-frequency noise and low-frequency baseline oscillations can be removed. Bioluminescence traces (Ts = 1. 67, N = 71) with increasing small molecule concentrations were fit with a damped sinusoid using the method described in a previous section. Because the small molecules were toxic at very high concentrations, experiments were removed from further analysis where R2 < 0. 80 (S1 Fig). A previously published mathematical model of circadian rhythms [28] was used to predict the effects on population damping rate from the dose-dependent small molecule experiments. The parameters used to capture the effects of each small molecule were the same as described previously [13]. The model was converted to a stochastic biochemical system and subsequently simulated using StochKit2 [29] (via GillesPy, https: //github. com/JohnAbel/gillespy). Population-level rhythms were found by taking the average of 1,000 noninteracting oscillators, starting at identical initial conditions. The only parameter left unspecified by the deterministic model was the cell volume, Ω, which controlled the amount of noise in the system. For each Ω, a damped sinusoid was fit to the population-averaged state trajectory. An R2 value was calculated for each fit, taking into account all eight state variables. We analyzed the data and annotations for the 111,743 wells (Ts = 2, N = 72) in the Zhang et al. , 2009 screen [12]. Fits for which the R2 < 0. 80 were discarded. The natural logarithm of the amplitude was used, since it more closely resembled a normal distribution and was on a similar scale to the damping rate. Plate to plate variation, as shown in S3 Fig, was more severe than variation on a well-to-well basis, S4 Fig. Parameters were therefore normalized on a plate-by-plate basis using a robust z-score [30]: z R, i = p i - M (p i) M (| p i - M (p i) |), i ∈ { 0, ⋯, P - 1 } where M (⋅) denotes the median of a vector, and pi contains all the points in one plate, and 𝒫 is the number of plates. We removed outlier points prior to calculating the moments of the distributions, Pearson correlation coefficients, and performing the multivariable linear regression. Outliers were defined as points that contained a z-metric (in either period, phase, amplitude, or damping rate) with an absolute value greater than eight. We chose the “control” wells to be those that contained no siRNA, as these proved to be more numerous than those containing reference siRNA perturbations and were clustered similarly to the highest-density regions of the perturbed fits. All computations were performed using Python. Code used to perform the analysis and produce the figures in this manuscript can be found online at https: //github. com/pstjohn/decay-methods.
  While both intrinsic and extrinsic noise sources can contribute to population-level damping, intrinsic noise is thought to play a more significant role in circadian systems [2]. We therefore first sought to determine whether changes to intrinsic stochastic noise alone are sufficient to explain population-level changes in damping rate. To do this, we calculated noise characteristics from experimental data on individual PER2: : LUC fibroblast cells [22]. Cells were sorted into two groups, a low-noise group and a high-noise group, based on the relative high-frequency noise, period variability, and amplitude variability present in each trace. Example rhythms from cells in both groups are shown in Fig 1A. Because the cells were not synchronized at the start of the recording, this effect is replicated in silico by shifting each series in time to align their start phases. Population-level bioluminescence traces were then found by averaging the cellular PER2: : LUC signal in each group. Both populations displayed averaged rhythms that resembled a damped sinusoid, similar to those seen in bioluminescence recordings of entire cell cultures. Fitting the averaged expression of each group with a damped sinusoid revealed that the low-noise group also had a lower damping rate (Fig 1B). The significance of this difference was confirmed via a bootstrap analysis (Fig 1C), where cells were randomly assigned in each bootstrap trial to either the low-noise or the high-noise group. Unlike from using single-cell imaging, inferring stochastic noise from the desynchronization rate of population-level recordings can be applied to existing and future high-throughput circadian screens. We demonstrate the insights that can be gained from such an approach by analyzing the publicly available genome-wide siRNA screen from Zhang et al. , 2009 [12]. The results of fitting a damped sinusoid to each of the 111,743 bioluminescence trajectories is shown in Fig 4, in which 86% of fits had an R2 > 0. 8. Since sinusoidal parameters can only be confidently inferred for fits in which the R2 was sufficiently high, wells were removed from further analysis if R2 < 0. 8. Additionally, of the fits with a high R2 value, only a small minority (0. 1%) had a negative damping rate. This trend supports the assumption that intercellular synchronization is unlikely in cultured U2OS cells. We next checked how parameters varied on a plate-to-plate and well-to-well basis. Well-to-well variation was relatively absent, aside from expected variation due to long- and short-period controls (S4 Fig). Fits were normalized to remove plate-to-plate variation (S3 Fig) using a robust z-score [30]. Additionally, we separated wells into a “perturbed” category and “control” category, depending on whether or not the well contained an siRNA perturbation. As we show in Fig 5, all fitted parameters displayed normal-like distributions, in which the control distributions showed much tighter clustering around the most likely values. Quantifications of the mean, variance, skewness, and kurtosis for each distribution are shown in Table 1. In the preceding sections, we have demonstrated that changes to single-cell intrinsic noise are sufficient to explain the observed changes in population-level damping. However, experimental work has shown that cell-autonomous fibroblast cells have a distribution of mean free-running periods [22]. Indeed, prior to the availability of single-cell data, studies explored the possibility that differences in mean periods served as the mechanism behind population-level damping [31]. While it is true that the dephasing of a group of oscillators can be caused by both variance in the mean period as well as cycle-to-cycle variability, intrinsic stochastic noise may play a more significant role. We show that there is greater variance in period on a cycle-to-cycle basis than on a cell-to-cell basis in cultured fibroblast cells (S5 Fig): individual cell period lengths have an average inner quartile range (IQR) of 3. 18 hrs, while cell-to-cell average period has an IQR of 1. 55 hrs. These results are independently confirmed by a previous study using Bayesian modeling, which found a standard deviation of periods within cells of 1. 43 hrs, and 0. 89 hrs across the population [32]. A similar result has also been observed in dispersed SCN neurons, suggesting that while both intrinsic and extrinsic period heterogeneity likely contribute to the dephasing kinetics, cell-to-cell differences are less severe than cell-autonomous noise [2]. It is also possible that damping rate changes due to siRNA or small molecule perturbation could be manifested through altering the system’s extrinsic noise. Such a change could be caused by an unequal uptake of siRNA or drug on a cell-to-cell basis, as has been demonstrated by a distribution of single-cell knockdown efficiency through flow cytometry [33]. This effect would increase cell heterogeneity and lead to faster dephasing kinetics. While differentiating between intrinsic and extrinsic noise sources from population-level data is possible in theory (S6 Fig), these differences are likely not identifiable from typical population-level data (S7 Fig). This distinction would likely be possible with single-cell level data, as has been done in other experimental systems [15]. However, such an experiment would likely not be amenable to high-throughput methods. Differences in damping rates are therefore best viewed as representative of changes to overall stochastic noise from both intrinsic and extrinsic factors. Since both types of noise are important to determining the overall function of population-level rhythms, damping rates are still a valuable method of quantifying stochastic noise.
  In this study we have shown that the damping rate of population-level circadian oscillations can be changed by genetic or pharmacological perturbations. As populations-level rhythms are determined by the coherence of many individual cells, desynchronization due to stochastic noise is a likely explanation for population-level damping. Using single-cell data, we showed that population-level damping rate is proportional to single-cell noise. Furthermore, we used a computational model to predict the changes in damping rate from two small molecules, demonstrating that changes to intrinsic stochastic noise are sufficient to explain the observed damping rate changes. We have described a method by which changes to stochastic noise can be estimated from population-level circadian bioluminescence recordings. An overview of the computational steps involved in our method are outlined in S8 Fig. While the method can be applied to existing experimental data, there are also practical considerations for the design of future experiments. Because the damping rate must be inferred from the time-varying amplitude, collecting bioluminescence data for longer time periods yields more accurate results and reduces the potential impact of initial transient regions. Additionally, a sampling rate that is high enough to confidently capture the peaks and troughs of gene expression is required—in this study, the 2 hr sampling window of the siRNA screen proved sufficient. While achieving such a rate is typically not difficult for bioluminescence experiments, it may limit the method’s applicability in experiments where samples need to be analyzed at each time point. We also note that there are many available tools for detrending and regressing time-series data. While we prioritized computationally efficient methods (which could be scaled to genome-wide screens), the best methods for any particular application may vary depending on the data. As high amplitude circadian oscillations are important for maintaining metabolic health, many studies have sought to find small molecule candidates that increase oscillatory amplitude [9]. In the search for circadian clock therapeutics, high-throughput methods are frequently used to screen for such drugs, often neglecting potential effects on cell-autonomous noise. While we have demonstrated a method by which the effects of small molecules on noise can be inferred from high-throughput methods, we have also shown that the potential improvement of clock robustness may be limited. Amplitudes of circadian rhythms may therefore be best increased by small molecule therapies that act transiently to synchronize peripheral oscillators. While such a method would require accurate alignment of drug administration to the correct circadian phase, a recent in silico study has demonstrated the potential effectiveness of such an approach in improving amplitudes in peripheral tissues [34]. Finally, the ability to extract an additional biologically relevant parameter from existing datasets will likely prove useful for many studies, as it allows further differentiation between perturbations that might otherwise have identical effects on bioluminescence rhythms.
summary: As most organisms exist in an environment that changes predictably with a 24-hour period, highly optimized genetic circuits turn on and off the production of key regulatory proteins to anticipate the day/night cycle. In humans, the demands of a modern society have required that we deviate from this evolutionarily prescribed sleep and feeding schedule, resulting in increased long-term risks of metabolic disease. There is therefore a desire to find pharmacological treatments that would restore the normal functioning of our circadian clock despite irregular behavioral schedules. One aspect of these treatments that is often overlooked in searching for candidate drugs is how these treatments might affect the accuracy of the circadian timing system. Recording the time of each cell is possible but difficult; as a result single-cell approaches cannot be scaled up to high-throughput searches. In this paper, we show that it is possible to estimate how much the noise of a system has changed by looking only at the averaged protein production of an entire population of cells. Such an approach allows us to analyze prior data from high-throughput screens, and show that the natural clock has been highly optimized to be both accurate and high amplitude.
title: Quantifying Stochastic Noise in Cultured Circadian Reporter Cells