Fruit fly activity was measured using the standard Drosophila Activity Monitor (DAM). Each fly is placed in an individual tube with food on one end and cotton on the other (Fig. 1a). An infrared beam crosses the tube in the middle in a perpendicular direction. When a fly interrupts the beam, the monitor receives a signal which is accumulated over time and sent to the computer every 20 s (Fig. 1b). Power spectra for activity data were calculated using maximum entropy spectral analysis (MESA) and Lomb–Scargle periodogram (LS) (Fig. 1c). Spectra obtained by each method show the expected 24 h peak. Additionally, there is a series of statistically significant peaks at smaller periods (\(T\)) with peak values from the two methods agreeing with each other to within 2 %. Unlike MESA, the Lomb–Scargle periodogram has an easily computable significance metric, which allows one to distinguish between significant and insignificant periodicities in the power spectra. Departing from the practice of filtering data prior to spectral analysis, we determine all power spectra directly from raw data. A number of past studies on ultradian rhythms adopted digital filtering as a standard procedure in their analyses [18, 21, 37–39]. However, as we demonstrate in this work (see “Methods”; Additional file 1: Fig. S9 and associated text), application of digital filters can irrevocably modify statistical properties of a time series and can give rise to artifacts in its power spectrum.
Figure 2a shows population averaged fly activity obtained from measurement of 14 flies in simulated light–dark conditions (LD) for 5 days. Drosophila daily activity has two distinguishable peaks. The morning peak (M) starts during the night and has a maximum when the light turns on and the evening peak (E) starts during the day and has a maximum when the light turns off. We constructed a model with a single fundamental period which reproduces these features of the activity. Our model consists of four normalized exponential terms with rates \(b_{MD} ,b_{MR} ,b_{ED} , \,{\text{and}} \,b_{ER}\), with subscripts denoting morning decay (MD), morning rise (MR), evening decay (ED) and evening rise (ER):
$$F(t) = \left\{ \begin{array}{*{20}ll} \frac{{e^{{b_{\small MD} T_{\small M} }} - e^{{b_{\small MD} t}} }}{{e^{{b_{\small MD} T_{\small M} }} - 1}} , &\quad 0 < t < T_{\small M} \\ \frac{{e^{{b_{\small MR} \left( {t - T_{\small M} } \right)}} - 1}}{{e^{{b_{\small MR} \left( {T_{\small 0} - T_{\small M} } \right)}} }} , &\quad T_{\small M} < t < T_{\small 0} \\ \frac{{1 - e^{{ - b_{\small ER} \left( {t - \frac{{T_{0} }}{2} - T_{\small E} } \right)}} }}{{1 - e^{{ - b_{\small ER} T_{\small E} }} }}, &\quad \frac{{T_{\small 0} }}{2} - T_{\small E} < t < \frac{{T_{\small 0} }}{2} \\ e^{{ - b_{\small ED} \left( {t - \frac{{T_{\small 0} }}{2}} \right)}} , &\quad \frac{{T_{\small 0} }}{2} < t < T_{\small 0} \end{array} \right.$$
In the proposed waveform, exponential terms with rates \(b_{MD}\) and \(b_{MR}\) form the M peak with width \(T_{M}\) while terms with rates \(b_{ER}\) and \(b_{ED}\) form the E peak with width \(T_{E}\) (Fig. 2a right). Together, these four exponents create a wave with \(T_{0} = 24\,{\text{h}}\). The Fourier transform of \(F(t)\) contains terms proportional to \({\text{cos}}\left( {2\pi \frac{{T_{0} }}{T}} \right)\) and \({\text{sin}}\left( {2\pi \frac{{T_{0} }}{T}} \right)\), which produce peaks in the power spectra at harmonics of the primary period \(T_{0}\) (see Additional file 1). To estimate peak heights, we analytically calculated Fourier transform only at harmonics \(T_{n} = T_{0} /n\) of the primary period. The analytic expression has multiple terms such as those containing \(\cos \left( {\frac{{2T_{M} n\pi }}{{T_{0} }}} \right)\), \(\sin \left( {\frac{{2T_{E} n\pi }}{{T_{0} }}} \right)\), \(e^{{T_{M} b_{MD} }}\), and \(e^{{T_{E} b_{ER} }}\) (see Additional file 1: equations 5–7).
Since a mathematical function that adequately reproduces fly locomotor data has been lacking, we tested different waveforms to simulate locomotion and interpret its power spectrum. In previous work, it was assumed that behavioral rhythms have the form of a square wave [37]. Although a square wave with period \(T_{0}\) produces multiple spectral peaks (Additional file 1: Fig. S1), it is able to interpret only odd-numbered harmonics in spectral data and does not completely mimic fly activity. We also tested a sawtooth wave, which is able to interpret most spectral peaks, but does not describe with high fidelity the shape of Drosophila locomotion data. Fruit fly circadian behavior is controlled by the activity of clock neurons, which often have exponential patterns of activation and deactivation [40]. The underlying exponential kinetics and the observed shape of activity data motivated us to build a model that consists of exponential terms.
In order to find similarities between the model and data and to look for the main periodicities in data, we first calculated the autocorrelation function. Autocorrelation with 1 min time lags was obtained by calculating the covariance of each time series with itself. The periods of the time series can be found from regularly appearing peaks of high correlation. Even though both the model and the data have a primary period of 24 h, the autocorrelations show strong 24 and 12 h periodicities (Fig. 2b). The 12 h periodicity results from the 12 h time interval between M and E peaks. Having obtained a basic level of similarity between data and model, we next analyzed them with the Lomb–Scargle periodogram to calculate a high-resolution power spectrum (Fig. 2c). Similar to the actual fly recordings, the proposed waveform with the primary period of 24 h shows multiple peaks in the power spectrum at values between 0 and 24 h with the 12 h peak dominating. The 12 h periodicity is significantly increased by the externally imposed 12-h light/12-h dark cycle. The model power spectrum not only reproduces prominent peaks from the data, but also the smaller peaks that appear near the prominent peaks. These small peaks, whose number and position can be predicted from the Dirichlet kernel, arise from the fact that we work with finite discrete data (Additional file 1: Fig. S8). For both the experiment and the simulation, the circadian peak is ~24 h.
We noticed that the majority of peaks that appear in the power spectra of wild-type flies measured in LD occur at multiples of the primary period \(T_{0}\). To test if the secondary peaks result from the externally imposed light/dark cycle, we measured fly (N = 29) locomotion in constant darkness (DD) for 5–7 days. This analysis revealed that together with the 24 h peak reflecting the circadian clock, the power spectrum still has the same additional peaks that appear at multiples of \(T_{0}\) (Fig. 3a, top), suggesting that the additional periodicities in the power spectra of fly activity are simply harmonics of the endogenous circadian period \(T_{0}\).
To further test our assumption that multiple spectral peaks of fly locomotion all result from the circadian clock, we used different circadian period mutants of Drosophila melanogaster. We analyzed per
S (N = 22) and per
L (N = 19) flies, with average circadian periods of 19 and 29.5 h (Additional file 1: Fig. S10), respectively, for the presence of secondary peaks in their power spectra. Activity of both genotypes was measured in constant darkness for 7 days and analyzed with Lomb–Scargle periodogram. Power spectra of the clock mutants show secondary peaks that, similar to wild-type, were found to be at multiples of the mutant \(T_{0}\) (Fig. 3a, middle and bottom panels). Almost all peaks in the power spectra of the different genotypes line up after rescaling the period axis with the circadian period. These results predict that elimination of the clock should abolish all rhythmicity in the power spectrum. We tested this prediction both by measuring circadian null mutants and by rendering wild-type flies arrhythmic by placing them in constant light (LL). Power spectra of per
0 (N = 38) and clk
Jrk (N = 23) clock mutants that do not show circadian rhythmicity in constant darkness were determined from locomotor recordings of individual animals (Fig. 3c; Additional file 1: Figs. S2, S4) [22, 41]. The majority of per
0 flies (N = 28) do not show statistically significant peaks at the p = 0.05 level between 2 and 35 h in the power spectra. Most clk
Jrk flies (N = 17) also do not have circadian rhythms and almost all peaks seen in the power spectra of wild-type flies between 2 and 35 h are absent in the power spectra of these mutant flies (Additional file 1: Fig. S2). We also measured activity of yw (N = 16) flies in LL and analyzed their power spectra (Fig. 3c; Additional file 1: Fig. S5). Most flies (N = 13) appear arrhythmic, showing no significant periodicities between 2 and 35 h. These results together provide strong support for the majority of spectral peaks resulting from the circadian clock producing a non-sinusoidal oscillation in fly behavior.
To test our model’s reliability, we compared peak positions in the power spectra for wild-type and shifted-period clock mutants of Drosophila to peak positions in the power spectra of model with different primary periods (Fig. 3b). For each genotype we constructed a model with \(T_{0}\) matching the circadian period in the data, yielding harmonics at values \(\frac{{T_{0} }}{1},\frac{{T_{0} }}{2}, \ldots ,\frac{{T_{0} }}{n}\) in the power spectra, which we then used to interpret peaks in the data. We used only peaks higher than the significance level of 0.005 in order to exclude Dirichlet kernel peaks from the analysis. For all three genotypes our simple model is able to identify more than 88 % of peaks in the data to within 10 % error (see “Methods” for details).
We next determined the exponents \(b_{MD}\) − \(b_{ED}\) for wild-type flies in DD and LD. In our model, these exponents define the shape of the M and E peaks. The model parameters were obtained from the power spectra of activity data (Fig. 4a). Spectra were fitted with an analytical expression \(H(T_{n} )\) obtained by calculating the square of the Fourier transform of \(F(t)\). The square of the Fourier transform yields peak heights \(H(T_{n} )\) at harmonics \(T_{n}\) of the primary period \((T_{n} = T_{0} /n)\) [see Additional file 1: equations (5)–(8)]. In DD data, \(T_{0}\) was determined from the peak at the circadian frequency. Since the fitting procedure is sensitive to the initial choice of parameters, as an initial guess we used parameters from a preliminary fitting of activity data with the model. The exponents obtained from the fits shown in Fig. 4 are \(b_{MD} = - 0.09\,{\text{h}}^{ - 1}\), \(b_{MR} = 0.34 \,{\text{h}}^{ - 1}\), \(b_{ER} = - 0.05\,{\text{h}}^{ - 1}\), and \(b_{ED} = 1.83\,{\text{h}}^{ - 1}\) for LD, and \(b_{MD} = 1.4 \,{\text{h}}^{ - 1}\), \(b_{MR} = 0.24\,{\text{h}}^{ - 1}\), \(b_{ER} = 0.50\,{\text{h}}^{ - 1}\), and \(b_{ED} = 6.2\,{\text{h}}^{ - 1}\) for DD. The analytical expression also produced good fits for the activity power spectra with the average peak height fit error of less than 10 % (Fig. 4b) (see “Methods” for details). Final values of the parameters determined from fitting the power spectrum were then used to construct a model for the activity recordings (Fig. 4c). For the measured wild-type flies, the constructed model shows good fit for locomotor activity with the average rate constant magnitude ~1.2 ± 2 h−1 (mean ± standard deviation). The parameters \(b_{MD}\) and \(b_{ER}\) obtained from fitting data can be either positive or negative, which means that the exponential terms with \(b_{MD}\) and \(b_{ER}\) in \(F(t)\) can be either concave or convex. Interestingly, we find that the parameters \(b_{MR}\) and \(b_{ED}\) are always greater than 0 and therefore corresponding terms are always concave. It should be noted here that our model does not impose any restrictions on the numerical values of the parameters.
While the predictive power of the model is reliable for \(T > 2 \,{\text{h}}\), for \(T < 2 \,{\text{h}}\) our model strongly deviates from data (Fig. 4a, inset). The calculated model peak heights \(H(T_{n} )\) predict that for period values less than \(2\pi /b\), where \(b\) is the largest exponent, typically \(b_{\text{MD}}\) or \(b_{\text{ER}}\), the peak heights are proportional to \(T^{4}\) (see Additional file 1). However, the data show weak sensitivity to T in this regime. There are a few factors that can affect the power spectrum of fruit fly activity at low period values. The first factor is that the locomotion measurement reports fly movement only when it crosses the middle of the tube. This lack of spatial resolution causes the system to miss small-scale movements which happen without the fly crossing the beam. A second factor is the Nyquist frequency. In any measurement the Nyquist frequency is equal to half of the sampling frequency and imposes a lower limit on the periods one can search for in a given time series. Thus, limitations imposed by both the measurement method and by the Nyquist frequency likely result in detectable discrepancies between the data and the model for low \(T\). For these reasons, all our model-based predictions are restricted to T > 2 h.
In order to understand what factors affect the model parameters, we determined their values for different clock mutants and tested for their relation to the circadian period \(T_{0}\). Wild-type, per
S, and per
L animals introduced above together with tim
UL flies (N = 11) with average T0 ~ 27 h were used in these analyses (Additional file 1: Fig. S10). Given our mathematical description of fly activity, we hypothesized two possibilities: one, in which the rate constants do not vary but \(T_{M}\) and \(T_{E}\) adjust with \(T_{0}\) and another, in which both sets of parameters change with \(T_{0}\). Interestingly, the data show that the rate constants \(b_{MD}\), \(b_{MR}\), \(b_{ER}\), \(b_{ED}\) do not depend strongly on \(T_{0}\) (\({\text{Adj }}R^{2} < 0.1\) for linear fits), whereas the parameters \(T_{M}\) and \(T_{E}\) that determine width of the morning and evening peaks, increase with \(T_{0}\) (Fig. 5a). Assuming the measured parameters represent characteristics of underlying biological processes, the robustness of the rate constants suggests that the pace of the clock likely does not alter the kinetics of these processes. On the contrary, tight association between \(T_{M}\) and \(T_{E}\) and clock speed suggests that the clock may regulate when these processes are initiated and terminated. The rate constants typically range between ±5 h−1 for \(b_{MD}\) and \(b_{ER}\), and between 0 and 1.5 h−1 for \(b_{MR}\) and \(b_{ED}\), with average \(\overline{{b_{MD} }} = 1.64 \pm 2.3 \,{\text{h}}^{ - 1}\), \(\overline{{b_{MR} }} = 0.78 \pm 0.77 {\text{h}}^{ - 1}\), \(\overline{{b_{ER} }} = 2.19 \pm 3.2 {\text{h}}^{ - 1}\), \(\overline{{b_{ED} }} = 1.2 \pm 1.6 {\text{h}}^{ - 1}\) (mean ± standard deviation). The mean ± standard deviation of the magnitude of the four rate constants from all flies tested in constant darkness is \(1.5 \pm 2.0 {\text{h}}^{ - 1}\). Lastly, linear fits of the \(T_{M\left( E \right)}\) versus \(T_{0}\) data reveal \(T_{M} (T_{0} ) \approx 0.19T_{0}\) and \(T_{E} (T_{0} ) \approx 0.29T_{0}\) (\({\text{Adj }}R^{2} = 0.90\) and \(0.95\), respectively) (Fig. 5a, bottom panels).
Based on the observed independence of the b parameters and the linear dependence of the \(T_{M(E)}\) parameters on the circadian period, we argued using the model that the average amplitude of daily activity must also increase with \(T_{0}\) (Fig. 5b). To test this prediction, we measured for each fly an activity amplitude, h, by averaging M and E peak heights in each recording. A plot of h vs. \(T_{0}\) shows that amplitude of activity indeed increases with lengthening of the circadian period (Fig. 5c). The form of our model function \(F(t)\) suggests that h and \(T_{0}\) can be related according to \({\text{h}} = c(1 - e^{{ - kT_{0} }} )\), which in the linear approximation becomes \({\text{h}} \approx CT_{0}\), with c, k and \(C\) as fit parameters. Due to the large scatter in our data, the simpler linear expression is statistically favored over the exponential function, yielding \({\text{h(}}T_{0} )\approx 1.4T_{0}\) (Adj \(R^{2} = 0.85\)). The surprisingly good agreement between the data and the prediction of our simple model further validates the mathematical description of fly locomotion pattern proposed in this work.