Circadian parameters mesor , circadian period , amplitude and acrophase for subjects A, B and C according to cosinor analysis, Fourier filter, discrete wavelet transform (DWT) using the Daubechies  mother wavelet and singular spectrum analysis (SSA) with parameter . Presented are average values (mean), standard deviation (SD) and the coefficient of determination .
In optimal conditions of youth and health, most—if not all—physiological systems obey regular circadian rhythms in response to the periodic day-night cycle and can be well described by standard techniques such as cosinor analysis. Adverse conditions can disturb the regularity and amplitude of circadian cycles, and, recently, there is interest in the field of chronobiology to quantify irregularities in the circadian rhythm as a means to track underlying pathologies. Alterations in physiological rhythms over a wide range of frequency scales may give additional information on health conditions but are often not considered in traditional analyses. Wavelets have been introduced to decompose physiological time series in components of different frequencies and can quantify irregular patterns, but the results may depend on the choice of the mother wavelet basis which is arbitrary. An alternative approach are recent data-adaptive time-series decomposition techniques, such as singular spectrum analysis (SSA), where the basis functions are generated by the data itself and are user-independent. In the present contribution, we compare wavelets and SSA analysis for the quantification of irregular rhythms at different frequency scales and discuss their respective advantages and disadvantages for application in chronobiology.
- singular spectrum analysis
- spectral analysis
- Fourier analysis
Circadian rhythms are physical, mental and behavioral variations that follow an approximately 24-hour cycle in response to the periodic alternation between day and night. In the last few decades, it has become well established that most—if not all—physiological systems exhibit these regular circadian rhythms, and it has been discovered that they are largely controlled by a central clock and several peripheral oscillators . More recently, it has been observed that adverse conditions, such as “healthy” and pathological aging, illness, stress and medication use, can disturb the regularity and amplitude of the circadian rhythm. Consequently, the focus in the field of chronobiology is shifting from a description of periodic cycles to the quantification of irregularities and the study of the mechanisms underpinning their disruption and normalization [2, 3, 4, 5].
The traditional method to analyze circadian rhythms is cosinor analysis, which quantifies the circadian 24 h cycle, and/or other specific infra- or ultradian periodic cycles, by means of examining the degree of “fit” between the experimental data and a user-defined model consisting of a superposition of cosine functions [6, 7] and allows to calculate the
Another disadvantage of the cosinor method is that it is unable to measure
The purpose of the present contribution is to illustrate how Fourier-based spectral analysis, wavelets and data-adaptive methods describe irregular circadian rhythms and rhythm fragmentation. Of the data-adaptive methods mentioned, in the present work, we prefer SSA because of its closeness to standard Fourier analysis and the availability of graphical tools such as the scree diagram that can be interpreted as a generalization of the well-known Fourier power spectrum. We will illustrate the advantages and disadvantages of these methods using selected 2-week continuous actigraphy time series for three nonsymptomatic control subjects of the public database of Ref.  (see Figure 1). These actigraphy time series show the number of movements per minute for a total duration of 20,160 minutes (h). The time series were chosen upon visual inspection: subject A is a young female adult (23 years old) with a regular circadian rhythm, subject B is a young male adult (22 years old) with an irregular circadian rhythm and subject C is an older male adult (82 years old) with a regular circadian rhythm.
2. Cosinor analysis
The traditional method to study the periodic aspects of circadian rhythms is cosinor analysis [6, 7]. The cosinor approach is based on regression techniques and is applicable to equidistant or non-equidistant time series of discrete data points:
The procedure consists of fitting a continuous periodic function to time series :
where is the average value or mesor around which the function oscillates, is the period, is the amplitude and is the phase which defines the value at which the function begins at the start of the monitoring at . By taking into account trigonometric rules, the cosinor function of Eq. (2) can be rewritten as , where and are amplitudes and where the phase is given implicitly by the superposition of the sine and cosine functions. Function is fitted to the data by minimizing the summed square residual errors for all data points , and the value of period for which the amplitude maximizes can be considered to be the circadian period (see Figure 2). For organisms exposed to the natural day-and-night cycle, this period can be expected to be h = 1440 min. Once chosen period , the other circadian parameters and are determined as well . Phase does not give any physiological information on the monitored individual, because time series may start at an arbitrary time of the day. Instead, a more interesting variable is the acrophase , which can be defined as the time of the day when the circadian cycle obtains its maximum, with respect to a fixed moment in time which is the same for all subjects, e.g. taking midnight as a reference, and which can be expressed as hours and minutes (hh:mm), or alternatively, as an angle (, taking into account the relation h), relative to this reference time. The fitted function may be generalized to include more than one period, , where the sum usually runs over a small number of different periods. In the present case, Figure 2 shows many ultradian (min) and infradian periods (min) that have nonzero amplitudes, but in the present case, there are no other periods than the circadian period h that are clearly distinguishable from the neighboring values to warrant their inclusion in the model function . An important result in cosinor analysis is the coefficient of determination , which compares the variance of the residual errors around the fitted model to the variance of the time series around its average value:
such that is a measure for the fraction of the variance of the time series that can be explained by model .
Results for the circadian parameters of the cosinor model function
3. Fourier spectral analysis
Fourier spectral analysis makes the supposition that the fluctuations of time series with may be interpreted as the superposition of independent harmonic oscillators, where is the Nyquist frequency and where each harmonic oscillator corresponds to a periodic function (see Ref. ):
which establishes a link between the time domain and the frequency domain . Here, the first term corresponds to or and is the direct current (DC) term around which the other terms oscillate and may be interpreted as the equivalent of the mesor of cosinor analysis. One of the most important results of Fourier spectral analysis is the power spectrum, which gives the power as a function of frequency . The variance of time series is given by
and Parseval’s theorem establishes that the variance in the time domain is identical to the variance of the oscillations of all components around the DC term in the frequency domain, i.e.
which allows us to interpret the power of the component with frequency as a
Results from the Fourier spectral analysis of subjects A, B and C are shown in Figure 3. The variance of the time series is Var(A) = 107,417, Var(B) = 110,405 and Var(C) = 28,725. Power spectra as a function of frequency in linear scale show for all subjects a dominant peak at the circadian frequency of oscillations during the 2-week duration of the time series. Apart from the circadian frequency, for subject C, only discrete peaks at higher harmonics are visible, whereas for subject A, many low-frequency components are present, which are even more predominant in the case of subject C. logarithmic scale allows to focus not on the dominant peaks but on the continuum of spectral contributions at a wide range of frequency scales, which for subjects A and B show an approximate scaling with spectral exponent and which for subject C appears to be truncated at about oscillations, below which the power spectrum flattens out. A drawback of the power spectrum representation as a function of is that it has a lot of dispersion because of which it may be difficult to make an accurate estimate of the value of the spectral exponent . The power spectrum can be re-ordered in the shape of a so-called
4. Wavelet analysis
4.1. Continuous wavelet transform (CWT)
The basic idea of wavelets is to decompose a time series in terms of a time-frequency set of orthonormal functions . The continuous wavelet transform (CWT), also called analytic wavelet transform (AWT), is a measure of similarity (in the sense of similar frequency content) between a basis wavelet function , called “mother wavelet”, and the signal itself. The resulting CWT depends on (see Figure 6 for a representation of some popular mother wavelets), but an appropriate choice of basis functions allows to analyze nonstationary time series. The Morlet mother wavelet has been suggested to be one of the most adequate bases for spectral studies [23, 24]. To evaluate CWT, the mother wavelet is translated to be centred at and scaled by a factor :
where the asterisk denotes the complex conjugate, and the result is the scalogram , which represents at each time the instantaneous period (represented by scale ) and the instantaneous intensity or power of the signal . At high frequencies, CWT has a good scale resolution but a poor frequency resolution, while at low frequencies, the frequency resolution is improved but time resolution is lost. One of the most important properties of a scalogram is the so-called ridge, which corresponds with the dominant behavior of the time series and in the present case can be expected to be the circadian rhythm. The extraction of the time-series component corresponding to the ridge can be not obvious, in particular when there are multiple ridges, and is the topic of the data-adaptive method of nonlinear mode decomposition (NMD) [23, 24]. The present discussion of CWT will be limited to the scalograms.
The left-hand panels of Figure 7 show the scalograms of the time series of subjects A, B and C using the Morlet mother wavelet. In the three cases, there is a high-intensity ridge near min, but the intensity and the instantaneous value of vary in time, reflecting the time and amplitude variability of the circadian rhythm. At ultradian frequencies, min, there are fringes of higher and lower intensities that alternate, reflecting day and night with high and low activities, respectively. The scalogram of subject C contains only a circadian ridge and fringes, whereas the scalogram of subject B shows a wide variety of competing and simultaneous features, and the scalogram of subject A is in between both the extremes. The right-hand panels show the average behavior of the scalograms over time, corresponding to a projection of the intensity values on the period axis, which are similar to the Fourier spectral analysis of Figure 4. For the three subjects, one can appreciate the dominant power of the circadian period of min and the scaling behavior at ultradian time scales. In absolute values, all scales have less power for subject C than for subjects A and B, but in relative values, the circadian cycle contributes much more for subject C than for the other subjects. Finally, to illustrate how CWT analysis depends on the wavelet basis, Figure 8 shows scalograms for the time series of subject A for other choices of the mother wavelet.
4.2. Discrete wavelet transform (DWT)
The discrete wavelet transform (DWT) starts with the partitioning of the signal into an approximation (smooth) and a detailed part, which both together yield the original signal itself. This subdivision is such that the approximation signal contains the low frequencies, whereas the detailed signal collects the remaining high frequencies. By repeatedly applying this subdivision rule to the approximation part, the details of increasingly finer resolution are then progressively separated out, while the approximation itself grows coarser and coarser. This procedure in effect offers a good time resolution at high frequencies and good frequency resolution at low frequencies, since it progressively halves the time resolution of the signal [14, 46]. The result of the DWT decomposition may depend on the wavelet basis chosen (see Figure 6 for an example of some popular mother wavelets), and the results may depend also on the number of scales into which the time series is decomposed. Ref.  suggests to use the Daubechies (4) mother wavelet for application of DWT to actigraphy data because of the discontinuous character of these time series.
Figure 9 shows the decomposition of the time series of subject A with DWT using a Daubechies (4) mother wavelet and a maximum number of scales, . The decomposition of the time series of the other subjects B and C is similar to the example shown for subject A. Visual inspection of the decomposition allows to identify the 10th scale, which oscillations about once every 24 h, as the circadian cycle. Figure 10 shows the circadian fit for subjects A, B and C. It can be appreciated that the fit reflects the variability in time and in amplitude of the time series. Corresponding circadian parameters are listed in Table 1. Figure 11 compares the partial variance carried by all 14 scales for the three subjects. It can be seen that both subject A and C present dominant circadian rhythms, whereas subject B is characterized by a large rhythm fragmentation at ultradian scales. Finally, Table 2 shows coefficients of determination for circadian fits for different choices of the DWT mother wavelet. Of all choices tried, Daubechies (4) appears to maximize the coefficient of determination .
5. Singular spectrum analysis (SSA)
SSA has been discussed in detail in a number of textbooks [16, 17, 18]; a short and very accessible introduction can be found in Ref. , whereas a larger and very complete review article is in Ref. . We have discussed the SSA method previously in Ref. . In brief, SSA can be explained as a three-step process: (i) the time series is transformed into a matrix which represents the underlying phase space of the time series, (ii) singular value decomposition (SVD) is applied to decompose this matrix as a sum of elementary matrices or—equivalently—to decompose the original phase space in a superposition of “subphase spaces” and (iii) each of the elementary matrices or “subphase spaces” is transformed back into a time-series component. Unlike Fourier analysis which expresses a time series as a sum of predefined sine and cosine functions, SSA can be considered to be
where is the time-series component, is the
Figure 12 shows some details of the decomposition of the time series of subject A. In the scree diagram, a trend mode , a mode which we will identify as the circadian (24 h) cycle and a mode which we will identify as an ultradian (12 h) cycle, can be distinguished with distinctive partial variances, and there is a long tail with components with similar partial variances that appear to obey a power law with . The weighted correlation matrix shows that the trend mode and the circadian mode are uncorrelated from the higher-order modes, but the -h mode does seem to have non-neglectable correlations with other components. The waveforms of before-mentioned modes confirm that they correspond to the trend, the circadian rhythm and a -h ultradian mode. Figure 13 shows the average correlation of the circadian mode with other SSA components, and broad minima can be observed for being equal or a multiple of the circadian periodicity, min. The more regular the circadian rhythm, the less correlated it is with the other time-series components and the better it can be isolated. Figure 14 shows the circadian fit with SSA for subjects A, B and C using min, and it can be appreciated that it describes the variability in time and amplitude of the time series. Corresponding circadian parameters are listed in Table 1. Figure 15 compares the scree diagrams for subjects A, B and C and is similar to Figures 4 and 7(d-e). In absolute values, it shows that time series of subjects A and B have comparable variances, which are much larger than the variance of the time series of subject C. In relative values, it shows that subject C has the strongest circadian rhythm and subject B the weakest. At the highest ultradian frequencies (), all subjects show a very similar behavior with a scaling. At the lowest ultradian frequencies (), scaling continues for subject A, whereas rhythm fragmentation is increased for subject B and decreased for subject C.
The interest of the field of chronobiology is shifting from a description of the periodicity of the circadian cycle to a quantification of deviations from regularity. The objective of the present contribution is to compare several methods in their description of irregular rhythms: the cosinor analysis, the Fourier filter, the continuous (CWT) and discrete wavelet transform (DWT) and the singular spectrum analysis (SSA). We are interested in irregular rhythms at the circadian time scale and rhythm fragmentation over a wide range of ultradian scales. Our aim is to illustrate the differences, similarities, advantages and disadvantages of the different methods using selected actigraphy time series.
We will first discuss the circadian time scale. According to the coefficient of determination
We will now discuss the ultradian rhythm fragmentation. The Fourier scree diagram, the CWT power spectrum and the SSA scree diagram suggest a trade-off effect. If subject A is taken as a reference, then subject C would seem to exhibit a very strong circadian rhythm associated with an important reduction of rhythm fragmentation at a wide range of ultradian scales, and, as a consequence, there is a flattening of the power law over the same frequency range, resulting in an overall rigid rhythm; subject B, on the other hand, is characterized by a reduced circadian contribution and an increased rhythm fragmentation over a wide range of ultradian scales, which leads to an increased power law slope over the same frequency range, resulting in an overall more random rhythm.
Of course, the three time series studied in the present contribution are not sufficient to draw any definite conclusions on the average values of the circadian parameters and their variability or on rhythm fragmentation at ultradian scales; therefore, a much larger statistical study is needed, but we have shown that circadian and ultradian scales can be studied within the same approach, and we hypothesize that partial variances are related over wide circadian and ultradian scales.
In recent years, there is a shift in interests in chronobiology where a larger emphasis is now put on an accurate quantification of irregularities of circadian rhythms and ultradian rhythm fragmentation to follow underlying pathologies. Wavelet analysis has probably been the method of choice to describe irregular rhythms at different time scales, but wavelet analysis has the drawback to depend on the choice of a mother wavelet which is arbitrary and user dependent. Data-adaptive time-series decomposition, where the basis functions are generated by the data itself, such as singular spectrum analysis (SSA) may offer an alternative. In the present contribution, we have shown that SSA is at least as versatile and accurate as wavelet analysis in the description and quantification of irregular rhythms at the circadian and ultradian time scales and may be a useful method to be adopted in the field of chronobiology.
Financial funding for this work was supplied by the Dirección General de Asuntos del Personal Académico (DGAPA) from the Universidad Nacional Autónoma de México (UNAM), grants PAPIIT IN106215 and IV100116 and in particular grant IA105017 in the context of which the Fourier filter was discussed into detail with Wady A. Rios-Herrera, Francisco F. de Miguel, Martha Yoko Takane Imay, David Serrano, Octavio B. Lecona and Silvia Diaz Gómez. We also gratefully acknowledge grants 2015-02-1093, 2016-01-2277 and CB-2011-01-167441 from the Consejo Nacional de Ciencia y Tecnologia (CONACyT). We are thankful for the Newton Advanced Fellowship awarded to RF by the Academy of Medical Sciences through the UK Governments Newton Fund programme. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript. Authors do not have any conflict of interest.