There is an array of numerical techniques available to estimate the period of circadian and other biological rhythms. Criteria for choosing a method include accuracy of period measurement, resolution of signal embedded in noise or of multiple periodicities, and sensitivity to the presence of weak rhythms and robustness in the presence of stochastic noise. Maximum Entropy Spectral Analysis (MESA) has proven itself excellent in all regards. The MESA algorithm fits an autoregressive model to the data and extracts the spectrum from its coefficients. Entropy in this context refers to “ignorance” of the data and since this is formally maximized, no unwarranted assumptions are made. Computationally, the coefficients are calculated efficiently by solution of the Yule-Walker equations in an iterative algorithm. MESA is compared here to other common techniques. It is normal to remove high frequency noise from time series using digital filters before analysis. The Butterworth filter is demonstrated here and a danger inherent in multiple filtering passes is discussed.

Physiological processes in almost all plants and animals have adapted to the cycles
in the environment, be they daily (circadian), tidal, lunar, synodic lunar monthly or
annual [

Analysis of time series may be simply done. In early work by Bünning on bean
plant leaf movements, periodicity was estimated by measuring peak-to-peak intervals
on chymographs that registered leaf position [

Circadian rhythms are studied in systems ranging from intracellular fluorescence
to complex behaviors such as running wheel activity; data acquisition and format
vary accordingly. For example, when studying the activity of an enzyme, the
variable may be continuous and an appropriate sampling interval must be chosen.
This must be rapid enough to avoid “aliasing” in the periodicity
region of interest. This occurs when the sampling interval is longer than the
period being recorded and is famously seen in old western movies when the spokes
of wagon wheels seem to be going backwards; sampling frequency must be no less
than twice the frequency of the cyclic process of interest. This constraint is the
Nyquist or foldover frequency [

Biological rhythm data are commonly not continuous and consist of unary events
unlike the record left on a kymograph by a bean plant. Here, other constraints
begin to play a role. Running wheel activity in mammals and the breaking of an
infrared light beam by

For the purposes of illustration of the techniques being discussed, we will
analyze a simulated data set. This is useful in this sort of discussion, as one
may know precisely what the parameters of the signal are. The series considered
here is a 23-h square wave with 20% white noise added consisting of 336 values
produced at half hour intervals for 7 days using our own software. I chose a
square wave as it is important to show that any analysis be effective for
waveforms deviating markedly from the badly overused sinusoid. A simple time plot
of the data is shown in Figure

Given a particular signal, even if it appears clearly rhythmic in a simple time
plot, it is important that an objective statistical test be employed to determine
if significant periodicity is present. Such a robust test is autocorrelation
analysis [

When computing the correlation coefficient, the output is normalized by dividing
by the variance of the complete data set, but this need not be so and the output
is then “covariance”, or the autocovariance function [

The autocorrelation function also yields a valuable way to quantify the regularity
of the signal both in terms of variation in period and the presence of noise. The
height of the third peak in this function, counting the peak at lag zero as one,
is taken as the Rhythmicity Index, or RI. This value relies on the decay of the
envelope of the function and is normally distributed so it may be used in
statistical analyses (review: [

It is useful to have a formal criterion for the significance of rhythmicity in
data. The 95% confidence interval for testing the significance of a given peak in
the autocorrelogram is 2/√N [

Beginning in the late 19^{th} century, the process of producing a spectrum
from digital time series was largely accomplished by Fourier analysis. Fourier
showed that any function showing certain minimal properties called the
“Dirichlet Conditions” can be approximated by a harmonic series of
orthogonal sine and cosine terms [

If our function consists of an ordered set of values x(

The exponential function consolidates the sine and cosine terms.

_{
ss
}

Fourier analysis has undergone considerable development and sees a great deal of
use in many fields, with chronobiology prominent among them; it has done yeoman
service. If the spectrum is calculated directly from data sampled at intervals, it
is termed the Discrete Fourier Transform or DFT. Fourier spectra are seldom
computed directly from the raw data however, rather they are produced from either
the autocovariance or autocorrelation functions. One argument for using the
autocovariance function is that the output is equivalent to partitioning the
variance in the signal by frequency and the area under the curve is the power
(review: [

Since the autcovariance and autocorrelation functions lose power with each pair of
points lost, usually no more than one third of the data are used to compute
correlation coefficients, adversely affecting the potential resolution in the
spectrum. To alleviate this, the rest of the function is padded out with zeroes.
This is an outright falsification of data points not in evidence, since there is
no reason to suspect these data points would all be zeroes. An added problem
occurs at the point where the zeroes abruptly start, since this abrupt
discontinuity will cause artifactual peaks in the spectrum called “side
lobes” owing to the Gibbs phenomenon [^{N} data points. Here again, the chances
of a data set containing an integer power of 2 points are slim, and again, zeroes
are added to pad the series out (Review [

Further data corruption can occur when tighter spacing of spectral estimates is
required. If the series is long, consisting of multiple cycles, this is usually
not a problem. However, when short data sets are at hand, as is commonly the case
with circadian rhythm work, there will be few cycles available. Fourier analysis
is based on harmonics and these are constrained. In practice, this means that the
spacing between estimates can be very wide. For a one-week long experiment with
data sampled at half hour intervals (as in the test data) and analyzed using a
simple DFT, spectral estimates are produced only for 22.4 and 24 h in the critical
interval between 22 and 25 h. This leaves enormous gaps with little chance that
the single estimate at 22.4 is even remotely close to the true period which is
23.0 h. Once again, it is straightforward to tighten up the interval between
estimates, but once again, zeroes are added with further problems in using false
data points for which there is no justification [

In the late 1960’s, John Parker Burg developed a new method for producing a
spectrum that tackles these problems [

The linchpin of this powerful technique is stochastic modeling. Time series evolve
in time according to probabilistic laws and there are a number of models that can
underlie such processes. One example is an autoregressive (AR) function (Review:
[

Where t is time, a is a coefficient derived from the data and Z_{t} is
white noise [

and, again, a, b, c,… are the model’s coefficients and p is the
order of the filter. These coefficients form the prediction error filter (PEF)
[

Where: _{k} is the set of PEF
coefficients.

The algorithm commonly used by us and others calculates the filter in an iterative
fashion and is based on the work of Anderson [

MESA has proven itself superior to ordinary Fourier analysis as it does not
produce artifacts from the various manipulations needed absent a model for the
function and both resolution and sidelobe suppression are superior to standard
Fourier analysis [

Biological signals can be notoriously non-stationary and noisy. This variation can
take the form of linear or nonlinear trends in amplitude, variations in period and
the waveform. As with any signal analysis system, MESA output can be improved by
conditioning the signal. It should be noted, however, that MESA is robust in the
face of such problems from the start. Incorporated into our MESA program is a
detrending step which fits a line by regression and subtracts it. This eliminates
linear trend and removes the mean. Mean removal is highly recommended, as this DC
component can obscure the rhythmic one if it is excessive [

We have had excellent success with Butterworth recursive filters [

Where X_{t} is the original data series and Y_{t} is the output
series. A and B are the filter coefficients: A = 9.656 and B = −3.4142. C
is the “gain” or amplitude change of the filter and equals 10.2426.
See [

MESA has seen notable success since first being implemented for use in biological
rhythms in the 1980s. MESA is useful for an extremely wide range of living
oscillatory processes. It was instrumental in discovering the presence of
ultradian rhythms in ^{0} and
^{-} mutations, which have no
overt circadian periodicity [

In summary, Maximum Entropy Spectral Analysis has proven itself to be a highly useful and versatile tool for the investigation of periodic biological phenomena.

A full explanation of the mathematics underlying MESA and the ways in which
algorithms have been implemented is beyond the scope of this paper. For those
wishing to explore these topics in detail, the author recommends the following:
For a good general introduction to the basic logic of MESA see Able’s
review [

All software used in this lab, including the FORTRAN source code, is available
free of charge from the author: dowse@maine.edu. A step by step annotated guide in
its use has been published by this author [

The author declares that he has no competing interests.