A Computational Method to Quantify Fly Circadian Activity
Summary
A method to quantify the main temporal features seen in fly circadian locomotor rhythms is presented. The quantification is achieved by fitting fly activity with a multi-parametric model waveform. The model parameters describe the shape and size of the morning and evening peaks of daily activity.
Abstract
In most animals and plants, circadian clocks orchestrate behavioral and molecular processes and synchronize them to the daily light-dark cycle. Fundamental mechanisms that underlie this temporal control are widely studied using the fruit fly Drosophila melanogaster as a model organism. In flies, the clock is typically studied by analyzing multiday locomotor recording. Such a recording shows a complex bimodal pattern with two peaks of activity: a morning peak that happens around dawn, and an evening peak that happens around dusk. These two peaks together form a waveform that is very different from sinusoidal oscillations observed in clock genes, suggesting that mechanisms in addition to the clock have profound effects in producing the observed patterns in behavioral data. Here we provide instructions on using a recently developed computational method that mathematically describes temporal patterns in fly activity. The method fits activity data with a model waveform that consists of four exponential terms and nine independent parameters that fully describe the shape and size of the morning and evening peaks of activity. The extracted parameters can help elucidate the kinetic mechanisms of substrates that underlie the commonly observed bimodal activity patterns in fly locomotor rhythms.
Introduction
The circadian clock is an endogenous biochemical oscillator with a period of approximately 24 hours and is almost ubiquitous in animals and plants1,2. The clock helps synchronize an organism's internal processes and behavior to the external light dark cycle. The genetic structure of the circadian clock has been widely studied since the 1960s using the fruit fly, D. melanogaster. In this insect, the core of the circadian clock consists of four proteins: PERIOD, TIMELESS, CLOCK, and CYCLE. These core components together with other molecules form a feedback loop that produces nearly sinusoidal oscillations of clock genes3,4. The circadian clock in flies is widely studied using multiday locomotor recordings where fly activity is detected with a single infrared beam crossing the middle of an individual tube5. A typical fly recording has a complex bimodal pattern with two well distinguishable peaks: Morning peak (M) that starts at the end of the night and has a maximum when lights turn on; and Evening peak (E) that starts at the end of the day and has a maximum when lights turn off6. Interestingly, the shape of such behavioral recording is very different from the simple sinusoidal oscillations observed at the molecular level, suggesting the action of additional mechanisms contributing to the observed temporal patterns. To better understand these hidden mechanisms, we have developed a computational tool that provides a quantitative description of the temporal patterns.
In our work, locomotor rhythms are defined in terms of a waveform that mimics fly activity pattern. Since simple sine waves cannot be used to model the observed rhythmic changes in activity, we tested various signal shapes to select the simplest one that captures all the salient features seen in the recordings. Fruit fly circadian behavior is controlled by the activity of clock neurons that often have exponential patterns of activation and deactivation7. The exponential dynamics and visual analysis of the data motivated us to build a model with exponential terms consisting of four exponents with nine independent parameters and closely resembling the fly activity pattern8. In addition to the locomotor data, we also analyze its power spectrum. Typical fly activity spectrum shows multiple peaks at harmonics T0/2, T0/3, etc., in addition to the expected fundamental peak at the circadian period T0. According to the Fourier theorem, only a pure sine wave produces a single peak in power spectra, while more complex waveforms show multiple spectral peaks at harmonics of the primary period (Figure 1). Therefore, given the non-sinusoidal temporal pattern in fly activity8, a multi-peak power spectrum of the data is mathematically expected and does not necessarily imply the presence of multiple periods of oscillation. Importantly, the power spectrum of the proposed model waveform also shows peaks at all harmonics of the primary period, similar to the fly locomotor recordings, thus underscoring the high fidelity with which our model describes fly data both in time and in frequency.
At time resolutions of a few minutes or less, fly activity data appears noisy, making it difficult to extract parameters directly from the raw data. Binning data into longer time intervals can decrease noise level, but, can alter the data in ways that can affect the estimation of model parameters. We therefore obtain the parameters from power spectra of the recordings, using an analytical expression for the expected power spectra calculated from the Fourier transform of the model function8 (see Additional File 1 of reference8). This approach of obtaining parameters from the power spectra yields accurate parameter values without any additional manipulations, such as binning or filtering, of the raw activity data. Mathematical details of the model and applications to wild-type and mutant data are described in reference8. The protocol presented here focuses on the step-by-step instructions to use the computational tool.
Protocol
1. Measuring Fly Locomotion Using Drosophila Activity Monitor (DAM)
NOTE: For more details see reference5.
- Prepare individual fly tubes with food on one end and cotton on the other. The end with food should be sealed to prevent the food from drying out.
- Put 5-6 g of fly food in a 50 mL beaker. Cut the food into small pieces so that it is easier to melt it.
- Connect 32 individual glass tubes with a rubber band.
- Melt the food in the beaker by heating it in a microwave oven for 10-15 s. Stop the microwave every 5 s and carefully shake the beaker to ensure equal melting of food.
- While the food is still liquid, put the prepared individual tubes in the beaker with food. Move the tubes up and down so they are equally filled.
- Allow the food to cool down and solidify for approximately 1 h.
- After the food is solid, remove the tubes with the food from the beaker.
- Seal the end containing the food using wax. First, carefully clean a tube using a paper towel, then press the tube against the wax. Visually check the quality of the seal, and if necessary, repeat the sealing again.
- Close the other end of the tube with cotton; cotton will permit air to go through while keeping a fly locked in the tube.
- Place a single fly in each individual tube and place the tubes in the DAM.
- Place monitors in an incubator that maintains constant temperature and humidity. Based on the experiment, set the proper light/dark conditions as follows.
- For light/dark experiments keep the flies in light/dark cycle for the whole experiment. Do not use the first day of measurements in the analysis.
- For constant darkness experiments, first keep the flies for two days in light/dark conditions for entrainment and synchronization of the clocks and then switch to constant darkness. Do not use the measurements from the first day of constant darkness in the analysis.
- Collect at least four days of data that can be used in the analysis.
NOTE: The DAM system will output a single file with a locomotion recording of all flies in the monitor.
2. Data Analysis
- Split the monitor output file into multiple single fly activity files; each file should be a single column '.txt' file with an individual fly locomotion measurement.
- Run 'ModelFitPS3.m' function in a Matlab command window with the following input parameters:
- For samplingrate, set the data sampling time interval in seconds. For example, if activity was measured every minute, enter 60 as the samplingrate.
- For bin_interval, set the time interval in minutes to which data will be binned for better visualization; the recommended bin interval is 20-30 min.
- For trend, enter "1" if the data show baseline trend and "0" otherwise; data with trend will be detrended first by fitting a second order polynomial and then subtracting it from the data.
- In the popup window select a single fly activity file.
Note: The first plot is data power spectrum, and not the familiar activity plot. From the plotted power spectrum, determine the primary period T0: Either click with the left mouse button on the circadian peak, or with right mouse button on the second harmonic peak (around T0/2). - On the opened data plot, check if the morning and evening peaks are well visualized. If not, change the bin_interval value by right clicking anywhere on the graph and inputting the new bin_interval value in the dialog box. The program will redraw the data with the new value of the interval. To accept the bin_interval value, left click anywhere on the graph.
- The program will redraw the data again and show the first five days of activity. On this plot, click on the first M peak that will be used in the analysis (sometimes it is necessary to skip one or two days).
NOTE: The program will redraw the graph starting from the picked morning peak. The blue and red lines will show the approximate position of the E peak and next day M peak, respectively, based on the period determined in step 2.4. - On the same graph, choose data for a preliminary fit of the data with the model: Click on the following points (in this order; note that the click location will be indicated with a red star on the bottom): (i) Top of M peak; (ii) End of the M peak; (iii) Start of the E peak; (iv) Top of the E peak; (v) End of the E peak; (vi) Top of the M peak of the next day.
- Note that the program now presents the power spectrum.
NOTE: The x axis is now given in frequency.- In the opened window with the power spectrum, pick points that will be used for fitting the analytical expression for the model power spectrum. The period detected in step 2.4 is marked with a red line. To pick fitting points, first roughly determine the primary period, similar to step 2.4. Then using the slider, fine tune the primary period value so that the fitting points (shown with red circles, will appear after moving the slider) are closed to peak values.
- After the visual peak selection, click "Accept" and the program will fit selected points with the analytical expression to calculate the model parameters.
- Note that the parameters and spectral fit error are saved to the file "model_fit_parameters.txt"; the program will additionally save 2 figures with fits of the locomotor data and its power spectrum.
Representative Results
The method presented here allows quantification of the main features in fly locomotion pattern. The quantification is achieved by fitting the activity data with a model that consists of four exponential terms:
The model has nine independent parameters that describe activity pattern. The parameters bMD,bMR,bED,and bER define the rates of morning decay (MD), morning rise (MR), evening decay (ED), and evening rise (ER), respectively. Parameter T0 defines circadian period, TM and TE define widths of M and E peaks, and HM and HE define heights of M and E peaks.
Values of the parameters are obtained from the unfiltered data in a few steps (Figure 2). First, a circadian period from the activity power spectrum (Protocol step 2.4) is determined. Then, preliminary parameters values are extracted from the locomotion recording by fitting it with the F(t) (Protocol steps 2.5-2.6). These values serve as initial guess for fitting the data power spectrum with an analytical expression derived in8 (Additional file 1 in reference8) by taking the Fourier transform of F(t):
where Tn = T0/n, with n = 1,2,3… and T0 is the circadian period (Protocol steps 2.7-2.8). By using the spectral fit, the model parameters are extracted from fly locomotion without filtering or binning. Fitting of the power spectrum produces final values for the model parameters, which are then used to construct the final form of F(t), i.e., the model of fly locomotor rhythms.
In the model, b parameters can be either negative or positive, which reflect the curvature of the corresponding exponential terms. The parameters bMD and bER are positive when the exponents are convex and negative when the exponents are concave, while the parameters bMR and bED are positive when the exponents are concave and negative when the exponents are convex. In general, for the exponents describing the M peak, the increase in the b parameters indicates slower rise or decay, and for the exponents describing the E peak, the increase in b parameters indicates faster rise or decay.
Figure 3 shows examples of fits obtained by this algorithm. Spectral peaks at T0 and harmonics from T0/2 to T0/10 are fitted (Figure 3A). Other harmonics typically have heights lower than significance level of 0.05 and are therefore not used in the analysis. Parameters obtained from the spectral fit are used to construct a model for fly locomotion (Figure 3B, red). The method works both for wild type flies (Figure 3, top panels) and for circadian mutants with altered period lengths (Figure 3, middle and bottom panels). For additional examples and connections between the model and underlying biology, refer to Lazopulo et al.8 Results from the fitting procedure are saved to a file in the following order: bMD,bMR,bER,bED, T0, TM/T0, TE/T0, HM, HE (Table 1). The program also saves to the output file the fitting error of the spectral fit, 'Err', which is calculated as a sum of squares of the residues (residual sum of squares) normalized by the square of the maximum of the data power spectrum.
Figure 1. Power spectra of different waveforms. According to the Fourier theorem, only a sinusoidal wave shows a single peak in a power spectrum, while more complex waveforms, such as square or sawtooth wave, show additional peaks at harmonics of the primary period. Please click here to view a larger version of this figure.
Figure 2. Schematic representation of the parameter-extracting algorithm. Final parameter values are obtained from fitting the power spectrum. Since the fitting procedure can be sensitive to starting values, the algorithm uses preliminary values calculated from the locomotion data and updates them using the power spectrum. The final parameter values are used to construct F(t) that best describes the data. Please click here to view a larger version of this figure.
Figure 3. Representative results of the method that quantifies fly locomotion using a model waveform. (A) Parameters of the model are extracted from fitting the first ten peaks of the power spectrum. (B) Obtained parameters are used to construct a model of activity. The method can be applied to wild type flies (top panels) or mutants with short (middle panels) and long (bottom panels) circadian rhythms. Please click here to view a larger version of this figure.
| fly genotype | bMD (1/min) | bMR (1/min) | bED (1/min) | bER (1/min) | T0 (min) | TM/T0 | TE/T0 | HM | HE |
| wild type | 0.0094 | 0.2077 | -0.0383 | 0.0606 | 1438 | 0.1528 | 0.0833 | 5.2238 | 8.7185 |
| short circadian mutant | 0.015 | 0.0086 | 0.5353 | 0.0227 | 1130 | 0.1404 | 0.2632 | 6.5481 | 7.3757 |
| long circadian mutant | 0.0069 | 0.0151 | 0.1035 | 0.9238 | 1701 | 0.2299 | 0.2644 | 7.2541 | 3.415 |
Table 1. Example of parameters extracted from fitting the power spectra shown in Figure 3. The program outputs the final parameter values to "model_parameters.txt" file as shown in the table.
Discussion
This work presents instructions for using a computational tool that provides a quantitative description of fly locomotion pattern. The tool fits locomotion data with a mathematical model consisting of four exponential terms that together describe the shape and size of the M and E peaks. The final values for the model parameters are obtained from fitting the power spectra of the data, where the use of the raw data can avoid artefactual effects that data binning or filtering can impose on parameter values. The model parameters can then be used to further study fly locomotion and the mechanisms that give rise to the activity peaks.
The quality of the fitting can be evaluated in two ways: by visual observation of the power spectrum fit and by using the fitting error (residual sum of squares) that is calculated together with the parameters. One of the key factors that influences the fit quality is the determination of the circadian period T0 (Protocol step 2.8). This value of T0 affects the estimation of peaks in the power spectrum and model parameters. Another key step that influences the fit quality is the selection of the starting point for a preliminary fit of the data (Protocol step 2.7), since TM and TE values are defined through this selection. The user can explore multiple fit results by choosing different days of activity for defining TM and TE.
Our model is focused on mimicking the well described clock controlled locomotor behavior with two prominent peaks of activity. This behavior is observed both in light/dark and in constant darkness conditions, which makes the approach applicable to these two types of experiments. However, the model does not take into account the light driven burst of activity observed in some light/dark experiments. Another limitation arises from the extraction of model parameters requiring the user input (Protocol steps 2.7-2.9). Since parameters TM, TE, and T0 are obtained from user input, each fly activity should be analyzed manually. Additionally, it is recommended to use the method with individual activities rather than with a population-averaged recording, since the averaging procedure hides many important individual differences in fly locomotion. Although it may seem cumbersome for large number of flies, by analyzing fly activities individually, the user obtains valuable statistical information about the degree of parameter variation within a group of animals. To overlook these benefits, it is possible to generate a population-averaged fly locomotion and analyze it to obtain average parameter values.
The model used in this work is an improved version of the waveform from our previous research8. Unlike the previous model, it has separate parameters for M and E peak heights, HM and HE, respectively. The new model predictably shows better agreement with fly activity data, since the morning and evening peaks typically have different heights.
This analysis is not limited to fly locomotion measured with infrared beam technique. Other methods, such as video-tracking, produce similar patterns of fly daily activity with M and E peaks showing exponential dynamics of rise and decay9. The tool discussed here can be readily applied to these alternative data sets as well.
The method discussed here quantifies fly locomotor recordings with the goal of connecting behavioral output to underlying mechanisms that regulate fly daily behavior and control the bimodal activity pattern. Several recent studies aimed at identifying neurons and substrates that shape M and E activity peaks10,11. The investigators manipulated the activity of certain neuronal groups and analyzed their role in the formation of the evening peak by observing changes in the height and phase of the E peak. In our model, parameters HE determines the height of the E peak, TE gives the peak width, and parameters bER and bED describe the peak shape. As outlined in reference8, in a biochemical context HE may indicate the measure of the levels of neuropeptides NPF and ITP, while bER and bED may represent the rates of their release and degradation. Together, these parameters can be utilized in future studies to better connect fly behavior with neuromodulatory signaling, towards an integrated description of the bimodal patterns in fly locomotor rhythms.
Disclosures
The authors have nothing to disclose.
Acknowledgements
We are grateful to Stanislav Lazopulo for help with the video content.
Materials
| Drosophila Activity Monitor | TriKinetics | DAM2, DAM5 | Measures fly locootion using single infrared beam |
| MatLab | Mathworks | Computing environment and programming language, MatLab should include Optimization and Symbolic Math toolboxes | |
| Drosophila melanogaster | per[S], per[L], iso31(wild type) | Our analysis can be performed with fly mutants of any circadian period |
References
- Pittendrigh, C. S. Circadian systems: general perspective. Biological Rhythms. II, 57-80 (1981).
- Zhang, E. E., Kay, S. A. Clocks not winding down: unravelling circadian networks. Nat Rev Mol Cell Biol. 11 (11), 764-776 (2010).
- Tataroglu, O., Emery, P. The molecular ticks of the Drosophila circadian clock. Curr Opin Insect Sci. 7, 51-57 (2015).
- Plautz, J. D., et al. Quantitative analysis of Drosophila period gene transcription in living animals. J Biol Rhythms. 12 (3), 204-217 (1997).
- Chiu, J. C., Low, K. H., Pike, D. H., Yildirim, E., Edery, I. Assaying locomotor activity to study circadian rhythms and sleep parameters in Drosophila. J Vis Exp. (43), e2157 (2010).
- Helfrich-Förster, C. Differential control of morning and evening components in the activity rhythm of Drosophila melanogaster–sex-specific differences suggest a different quality of activity. J Biol Rhythms. 15 (2), 135-154 (2000).
- Dautzenberg, F. M., Neysari, S. Irreversible binding kinetics of neuropeptide Y ligands to Y2 but not to Y1 and Y5 receptors. Pharmacology. 75 (1), 21-29 (2005).
- Lazopulo, A., Syed, S. A mathematical model provides mechanistic links to temporal patterns in Drosophila daily activity. BMC Neuroscience. 17 (1), 14 (2016).
- Donelson, N., Kim, E. Z., Slawson, J. B., Vecsey, C. G., Huber, R., Griffith, L. C. High-resolution positional tracking for long-term analysis of Drosophila sleep and locomotion using the "tracker" program. PloS ONE. 7 (5), e37250 (2012).
- Schlichting, M., et al. A Neural Network Underlying Circadian Entrainment and Photoperiodic Adjustment of Sleep and Activity in Drosophila. J Neurosci. 36 (35), 9084-9096 (2016).
- Guo, F., et al. Circadian neuron feedback controls the Drosophila sleep-activity profile. Nature. 536 (7616), 292-297 (2016).
