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 plants^1,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 genes^3,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 tube⁵. 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 off⁶. 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 deactivation⁷. 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 pattern⁸. In addition to the locomotor data, we also analyze its power spectrum. Typical fly activity spectrum shows multiple peaks at harmonics T₀/2, T₀/3, etc., in addition to the expected fundamental peak at the circadian period T₀. 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 activity⁸, 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 function⁸ (see Additional File 1 of reference⁸). 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 reference⁸. 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 reference⁵.

1. 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.
   1. 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.
   2. Connect 32 individual glass tubes with a rubber band.
   3. 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.
   4. 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.
   5. Allow the food to cool down and solidify for approximately 1 h.
   6. After the food is solid, remove the tubes with the food from the beaker.
   7. 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.
   8. Close the other end of the tube with cotton; cotton will permit air to go through while keeping a fly locked in the tube.
2. Place a single fly in each individual tube and place the tubes in the DAM.
3. Place monitors in an incubator that maintains constant temperature and humidity. Based on the experiment, set the proper light/dark conditions as follows.
   1. 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.
   2. 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.
4. 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

1. 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.
2. Run 'ModelFitPS3.m' function in a Matlab command window with the following input parameters:
   1. For samplingrate, set the data sampling time interval in seconds. For example, if activity was measured every minute, enter 60 as the samplingrate.
   2. 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.
   3. 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.
3. 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 T₀: Either click with the left mouse button on the circadian peak, or with right mouse button on the second harmonic peak (around T₀/2).
4. 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.
5. 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.
6. 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.
7. Note that the program now presents the power spectrum.
   NOTE: The x axis is now given in frequency.
   1. 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.
8. After the visual peak selection, click "Accept" and the program will fit selected points with the analytical expression to calculate the model parameters.
9. 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 b_MD,_bMR,_bED,_{and bER} define the rates of morning decay (MD), morning rise (MR), evening decay (ED), and evening rise (ER), respectively. Parameter T₀ defines circadian period, T_M and T_E define widths of M and E peaks, and H_M and H_E 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 in⁸ (Additional file 1 in reference⁸) by taking the Fourier transform of F(t):

where T_n = T₀/n, with n = 1,2,3… and T₀ 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 b_MD and b_ER are positive when the exponents are convex and negative when the exponents are concave, while the parameters b_MR and b_ED 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 T₀ and harmonics from T₀/2 to T₀/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.⁸ Results from the fitting procedure are saved to a file in the following order: b_MD,_bMR,_bER,_bED, T₀, T_M/T₀, T_E/T₀, H_M, H_E (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 T₀ (Protocol step 2.8). This value of T₀ 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 T_M and T_E values are defined through this selection. The user can explore multiple fit results by choosing different days of activity for defining T_M and T_E.

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 T_M, T_E, and T₀ 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 research⁸. Unlike the previous model, it has separate parameters for M and E peak heights, H_M and H_E, 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 decay⁹. 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 peaks^10,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 H_E determines the height of the E peak, T_E gives the peak width, and parameters b_ER and b_ED describe the peak shape. As outlined in reference⁸, in a biochemical context H_E may indicate the measure of the levels of neuropeptides NPF and ITP, while b_ER and b_ED 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.

Materials

Name	Company	Catalog Number	Comments
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