A mathematical model provides mechanistic links to temporal patterns in Drosophila daily activity

Biological oscillators with periods varying from seconds to years play important roles for most living organisms [1]. These oscillators have been divided into three groups by period length: (1) circadian oscillators, with periods close to 1 day; (2) ultradian oscillators, with periods less than 24 h; and (3) infradian oscillators, slower than circadian oscillators, with periods from a few days to a few seasons. Examples of ultradian rhythms can be found in the oscillation of body temperature in golden hamsters and in electroencephalogram measurements from the human brain [2, 3]. Infradian oscillations can be seen in animals, such as birds, that have annual migration cycles [4]. Behavioral oscillators form the group of oscillators that are responsible for adaptation to external environments, control of daily cycles of activity, annual reproductive rhythms, and migrations [5]. The most prominent and well-studied of the behavioral oscillators is the circadian clock [6–8].

The circadian clock is an endogenous biological oscillator with a period of approximately 24 h that controls daily activity and is found in most plants and animals. This clock helps to synchronize an organism’s internal processes and behavioral outputs with the daily light–dark cycle. The fruit fly Drosophila melanogaster is one of the primary model organisms in the study of behavioral rhythms and previous studies in this invertebrate have revealed the core components of the circadian clock that have been subsequently found to be conserved in other organisms, including mammals [9].

How these overt behavioral rhythms are generated has been an active topic of research for decades. One of the first theoretical models for the circadian system was proposed by Pittendrigh, who suggested that the circadian clock consists of multiple oscillatory components that maintain identical frequency and appropriate phasing through external entrainment [10]. An alternative mechanism was proposed by Pavlidis who relied on contemporary research showing that most known biochemical oscillators had periods which did not exceed a few minutes [11]. Pavlidis showed mathematically that a system of strongly coupled oscillators with short periods can produce a robust circadian rhythm and simulate many aspects of circadian behavior. The Pavlidis mechanism was supported subsequently by various mathematical models in which circadian rhythms result from a coupling of ultradian oscillators with periods of a few hours [12, 13]. In yet another theoretical study circadian rhythms were generated by coupling oscillators with periods of only a few seconds [14].

Behavioral studies of various animals claim to reveal multiple ultradian rhythms in daily activity [15–21]. These studies are commonly performed using spectral analysis of activity recordings. Power spectra typically show multiple peaks at periods of less than 24 h, which some studies have attributed to ultradian oscillators. In these behavioral studies short ultradian rhythms associated with foraging activity were observed in Siberian hamsters [15] and common voles [16]. The animals were studied under conditions disruptive to circadian rhythms, leading the authors to suggest that these ultradian rhythms are independent of the circadian clock. Other research has shown periods between 4 and 12 h in the activity of inbred strains of mice [17]; these rhythms were observed both in 12 h light/12 h dark conditions and in constant darkness. A more recent report from Blum et al. on dopamine-dependent oscillations in mouse activity with a ~4 h period has reinvigorated the discussion on ultradian rhythms in mammalian behavior [18]. In studies of fruit fly activity, ultradian rhythms were first claimed in per ⁰ flies measured under constant conditions [19, 20]. Circadian rhythms are abolished in per ⁰ animals due to a mutation in the period (per) gene [22]. Subsequently, ultradian rhythms were also reported in wild type and other clock mutants of Drosophila melanogaster [21]. Based on these data, the authors proposed that in fruit flies per couples ultradian oscillators and ultimately produces a circadian rhythm. However, even though it was not explicitly mentioned, power spectra in these studies were interpreted under the assumption that only one spectral peak results from the circadian clock and therefore the activity data have a sinusoidal shape. In reality, rhythmic behavioral data often have a non-sinusoidal shape in time.

The core of the Drosophila circadian clock consists of four proteins, PERIOD (PER), TIMELESS (TIM), CLOCK (CLK) and CYCLE (CYC), that form a negative feedback loop [7]. Although the molecular components differ, circadian clocks in mammals and plants also employ similar genetic architecture and, in each case, this conserved architecture produces clock gene oscillations that are sinusoidal in shape with period of approximately 24 h [8, 23–26]. These molecular oscillations converge on downstream behavioral circuits such as those for sleep, feeding and mating, and ultimately shape daily activity of organisms [3, 27, 28]. In Drosophila, temporal patterns in solitary activity are generally studied using locomotion measurements by detecting a fly crossing an infra-red light beam in the middle of a tube. A typical recording shows a non-sinusoidal time series with two distinguishable peaks: the morning peak (M), which starts during the end of the night time and has its maximum when the light turns on; and the evening peak (E), which starts during the end of the day and reaches its maximum when the light turns off. The M and E activity peaks are thought to be produced by two groups of pacemaker neurons, the M and E oscillator cells, and modulated by neuropeptides secreted by the pacemakers [29–32]. These two groups of neurons produce different neuropeptides. The neuropeptide produced in the morning oscillator cells is pigment dispersing factor (PDF) which has been shown to synchronize clock neurons and promote morning activity in fly behavior [33]. One of the neuropeptides produced by the evening oscillator cells is called the ion transport peptide (ITP) which is thought to influence fly activity around dusk [34]. A second possible source for the non-sinusoidal shape is the fly sleep homeostat, a feedback system that keeps track of sleep need in the animal. In mammals, the homeostat output has been modeled as repeating patterns of exponential rise and decay [35]. Assuming the fly sleep homeostat produces similar patterns, interaction of its output signal with that of the clock could shape the final behavioral output as a non-sinusoidal waveform, as seen in locomotor recordings. Indeed, that sleep–wake signals influence locomotor output is a standard assumption in fly studies which now routinely use locomotion as a read-out of sleep and wakefulness [36]. These possibilities suggest that understanding the temporal shape of the recordings could provide systems-level access to the biochemical pathways that influence fly locomotion patterns.

A mathematical description of the widely studied daily locomotion patterns in flies is currently missing and its absence has limited our interpretation of the data in both frequency and time domains. In frequency space, spectral data of locomotion are typically interpreted without considering the shape of the activity time series, although it is well-known from Fourier’s theorem that the form of the time series is a critical determinant of its power spectrum. For example, non-sinusoidal signals with a single periodicity, such as square or triangular waves, produce power spectra with multiple peaks, which correspond to the harmonics of the primary period. In the time domain, beyond a qualitative description of how M and E peak heights vary with simple perturbations, no biophysical description exists for the peculiar shape of fly locomotor activity. For instance, we do not know what biological processes control the rates of increase and decrease of activity around M and E peaks.

We address these limitations by constructing a simple waveform with a single period that closely resembles the shape of fruit fly activity. Our model consists of four exponential terms that generate a pattern resembling the M and E peaks. The power spectrum of the proposed waveform has peaks at harmonics of the primary period \(T_{0}\) and we demonstrate that the predicted spectral peaks can account for the multiple periodicities seen in the power spectrum of Drosophila locomotion. The data show that lesions in circadian genes that shift the circadian peak also shift the ultradian peaks, in accordance with our model. Genetic or environmental ablation of the circadian clock expectedly eliminates the circadian spectral peak but also reduces to noise levels the height of peaks in the ultradian range. Thus we weaken the competing view that circadian rhythms may arise from the coupling of ultradian oscillators by showing that the secondary spectral peaks in Drosophila locomotion, rather than being produced by independent oscillators, simply result from the non-sinusoidal shape of fruit fly activity. In the time domain, we determine the model exponential rate constants for wild-type and clock mutants of Drosophila. These rates which vary widely in magnitude between 0.024 and 14.5 h⁻¹, with mean value of 1.49 h⁻¹, are found to be independent of the pace of the circadian clock. Guided by these results, we propose that the exponents might come from the accumulation and release rates of neuropeptides important in regulating activity and sleep in the fly brain. Finally, we use the model to make quantitative predictions about candidate biological processes in general that could give rise to the rate constants and suggest future experiments to help identify the associated molecular substrates.

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 ⁰ (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 ⁰ 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⁻¹ (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 T₀ ~ 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⁻¹ for \(b_{MD}\) and \(b_{ER}\), and between 0 and 1.5 h⁻¹ 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.

While a large number of studies have exploited rhythmic patterns in fruit fly activity to query the status of the circadian clock, the broader question of what principal endogenous systems other than the clock shape these patterns, remain poorly understood. Here we attempt to fill this gap by formulating a simple mathematical model that accommodates the diverse locomotion patterns exhibited by Drosophila over the time-scale of multiple days. The model has two major components: an oscillatory component characterized by a single period and an exponential component defined by four rate constants. We show that these components are together sufficient to explain prominent features in both the power spectra and the time series of fly locomotor recordings. Using the model, we show that the position and height of over 88 % of locomotor spectral peaks can be accounted for as harmonics of circadian rhythms. We use this finding as evidence to support the view that in flies daily locomotor rhythms receive negligible input from endogenous ultradian clocks. Using the proposed model, we additionally extract values of the exponential rate constants in a variety of flies. We propose that these rate constants could represent kinetics of neurobiological processes that, such as the secretion of neuropeptides, have been implicated in modulating behavior. We suggest our data and model-based predictions open up fly behavioral recordings to more rigorous analysis.

Ultradian and circadian clocks are known to co-exist in animals but the extent of interactions among the two classes of endogenous oscillators is still unresolved. Previous work on this issue outlines three possible relationships between circadian and ultradian rhythms: i) a circadian rhythm is produced by the coupling of ultradian oscillators, and ultradian peaks seen in power spectra are components that break down from the circadian clock [19–21, 42]; ii) ultradian and circadian oscillators exist independently in the organism [15, 21]; and iii) ultradian rhythms come from the desynchronization of a population of circadian oscillators [21]. Using our model we can address the first two hypotheses. We have shown that multiple periodicities in power spectra result from a non-sinusoidal shape of Drosophila locomotion and appear at harmonics of the circadian rhythm (Fig. 3a, b). Since nearly all peaks in our fly power spectra come from circadian rhythms (Fig. 3c; Additional file 1: Figs. S2–S5), we suggest the first hypothesis, which states that circadian rhythms result from the coupling of ultradian oscillators, is unlikely. The assumption of the second hypothesis, that circadian and ultradian oscillators coexist, is possible; however we have shown that more than 88 % of peaks that stably appear between 2 and 35 h in the power spectra of fly locomotion result from the circadian clock. Our attempts to detect intermittent ultradian activity revealed such behavior in <5 % of the total recording time (Additional file 1: Fig. S7). These results therefore indicate that there are no robust ultradian oscillators with periods in this range detectable in fly activity. Our data are at odds with claims made by several previous examinations of locomotion [17, 19–21, 43]. The previous studies missed the harmonics likely due to a combination of disregarding the non-sinusoidal pattern of the raw time series, the inappropriate use of digital filters and the lack of proper statistical metrics. Our data, together with bona fide ultradian periodicities discussed in other behavioral contexts such as male courtship song [44, 45] or calcium oscillations in neurons [46], suggest that second hypothesis is most likely, and that there is minimal cross-talk between the two types of oscillators, at least at the level of behavioral output. In brief, our results contradict the previously espoused role of multiple ultradian oscillators shaping locomotion. Instead, the data support a parsimonious model in which a single circadian oscillator regulates all rhythmic modulations of fly locomotion.

Data on ultradian rhythms reported in locomotion of mammals can also be considered in relation to the circadian clock. Studies that addressed this question in Siberian hamsters and common voles concluded that the ultradian rhythms are generated independently of the circadian clock [15, 16]. Autonomy of mammalian ultradian rhythms was further affirmed by a recent work in mice that identified dopamine as a major regulator of and, possibly, the source of these ~4 h oscillations [18]. These previous data from mammals raises the question of why Drosophila locomotion does not exhibit robust ultradian oscillations. Our model, which incorporates only circadian rhythms, reliably predicts spectral peak heights and positions in rhythmic flies of different genotypes and genetic backgrounds and shows an absence of rhythms in circadian arrhythmic animals. Compared to mammals, the fly data suggest that if there is an ultradian signal that feeds into locomotion, it is either much weaker or has a period ≪2 h where our model does not strictly apply.

Beyond addressing the rhythmic component of locomotion, our model also captures the general temporal shape of fly locomotion. Despite widespread use of these time series in fly circadian studies, to our knowledge, proposed \(F(t)\) is the first description of the multi-day data in mathematical terms. In our mathematical model, the various shapes of these time series are quantified in terms of four exponential rate constants, \(b_{MD} - b_{ED}\), that control the slopes around the morning and evening activity peaks and the parameters \(T_{M} - T_{E}\), that control the widths of these peaks. Although one may assume that all model parameters should scale with the circadian period, our measurements reveal that while \(T_{M}\) and \(T_{E}\) increase with the circadian period, the exponents \(b_{MD} ,\,b_{MR}\),\(b_{ER}\), and \(b_{ED}\) intriguingly do not (Fig. 5a). Independence of the rates from the circadian period indicates that mutations in per and tim do not affect the rate constants \(b_{MD} - b_{ED}\), suggesting that the rate constants possibly represent clock-autonomous biological processes. This observation together with available biochemical data, leads us to hypothesize that the rate constants could represent processes like accumulation, release or degradation of neuropeptides, particularly those that are involved in the clock neuronal network. Neuropeptides pigment dispersing factor (PDF), ion transport peptide (ITP) and neuropeptide F (NPF) modify M and E clock neurons and influence locomotor patterns [34, 47, 48]. The abundance of these neuropeptides, in turn, are modulated by the clock [34, 49]. But their rates of accumulation in neurons, subsequent secretion, and degradation are presumably controlled by cell-biological machinery whose rate-limiting steps are set independently of the circadian program [50]. Therefore, we propose the \(b\) rate constants could represent average measures of some of these kinetic steps. Specifically, since PDF controls morning anticipation, which in our model is shaped by \(b_{MD}\) and \(b_{MR}\), we suggest these parameters may reflect kinetic steps related to PDF [33]. Similarly, because ITP reduces nocturnal activity and enhances diurnal evening activity, its kinetics may contribute to exponents \(b_{ER}\) and \(b_{ED}\) which form the E peak in our model [34]. Regarding temporal modulation of the levels of the neuropeptides, our model suggests that it likely results from the clock periodically initiating and terminating the key kinetic steps. In circadian rhythms, neuropeptide activity has so far been accessible largely through biochemical assays. Our analysis using \(F(t)\) may have established a potential link to the relevant neuropeptide kinetics through the behavioral rhythms that these neuropeptides are thought to modulate.

Alternatively, the model exponents could represent signals that feed into locomotor circuits from other endogenous systems such as those controlling fly sleep–wake states. Daily temporal patterns in sleep and activity in flies, like in mammals, are coordinated with major inputs from the circadian clock and the sleep homeostat [36]. Studies on human sleep patterns, in particular, have led to a model in which the sleep homeostat output is described as an exponentially varying signal [35]. In this model, the varying signal is the sleep electroencephalogram (EEG) power density that exponentially increases during daytime wakefulness and diminishes during nighttime sleep, producing a singular peak in the signal in the early evening. The average rate constants associated with the increase and decrease of EEG power are ~0.06 and ~0.24 h⁻¹, respectively. Interestingly, these are well within the range of values we find for the \(b\) parameters. Although output of the fly sleep homeostat has yet to be defined, it is reasonable to assume that its architecture may resemble that of the human homeostat given the substantial mechanistic overlap already uncovered between fly and human sleep [36]. Together, these similarities are consistent with a hypothesis in which the rate constants extracted from fly locomotion could be representative of the fly sleep homeostat.

Our mathematical model also permits a number of general predictions to be made about candidate neurobiological processes that may underlie our model parameters (Fig. 6). First, we propose that the candidate process can be described, for simplicity, in terms of linear or exponential growth-decay kinetics. Second, the processes are presumably sporadic when averaged over many neurons. Their random nature would imply they have noisy power spectra (Fig. 6a, middle row) and considering the processes are key modulators of behavior, their average stochastic nature should result in arrhythmic behavior in the absence of the circadian clock (Fig. 3c). Third, we suggest that the underlying processes have kinetics on the time-scale of tens of minutes. This constraint is based on the magnitude of the model rate constants, which are ~1.5 h⁻¹ on average, suggesting that the processes they represent persist for ~40 min. At the molecular level, this is a long time-scale and suggests that neurotransmitter activity or electrical pulses, which also have exponential kinetics but act on the time-scale of seconds or less, are unlikely to be direct contributors to the exponents in \(F(t)\). However, neuropeptide half-life in the brain and accumulation in synaptic boutons happen over span of ~30 min, making these substrates attractive candidates in our model [50]. Lastly, we suggest that on the time-scale of hours the circadian clock is likely the sole pacemaker (Fig. 6a, top row) to temporally gate the relevant neurobiological processes. The circadian gating imposes rhythm on otherwise stochastic processes and when integrated with additional downstream signals, ultimately produces rhythmic patterns in locomotion (Fig. 6a, bottom row). Related to the circadian gating are two additional features that the hypothetical processes might display (Fig. 6b). We suggest that changes in the period of the clock may cause not only parallel changes in the rhythm of the processes but, importantly, also changes in the amplitude of their oscillation. For instance, if \(b_{MD}\) and \(b_{MR}\) are indeed related to PDF activity, our model predicts that compared to per ^L animals, in per ^S animals both the period and the peak–peak amplitude of PDF oscillation is smaller. This prediction is based on our assumption that such changes in the amplitude of the underlying processes likely causes the correlation observed in our data between the average locomotor amplitude and the circadian period, which is predicted surprisingly well by our model (Fig. 4b, c). Together, these characterizations should facilitate experimental identification of the key neurobiological processes that our model parameters represent.

In summary, this work initiates a more mathematical dissection than is currently available of the temporal shape of fly locomotor recordings. We point out that such analysis can lead to experimentally verifiable predictions of what major biological processes regulate daily manifestations of behavior. As examples of such candidates, we propose accumulation and degradation of neuropeptides implicated in controlling sleep-activity in flies and suggest these kinetic steps are gated principally by the circadian clock. Critically, our mathematical model now opens up the recordings to quantitative scrutiny towards clearer mechanistic understanding of daily rhythmic behaviors in Drosophila and other organisms.