Quantification of Irregular Rhythms in Chronobiology: A Time- Series Perspective

Abstract

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.

1. Introduction

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 [1]. 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 circadian parameters of mesor, amplitude, period and acrophase [8]. However, data where the patterns are irregular, or where the statistical properties vary over time (nonstationary time series), such as having a dominant trend [9, 10, 11], or time-varying amplitudes, frequencies or phases [12, 13, 14, 15], are much harder or impossible to describe using models based on these periodic functions. Recently, more specialized techniques have been developed to study circadian rhythms; in particular, wavelets have been applied to study the irregular aspects of circadian rhythms [13, 14, 15]. Wavelets however are, as with cosinor analysis, model-based in the sense that the results obtained may depend on the particular wavelet basis functions (mother wavelet) selected by the user. Between the most recent developments in the field of time-series analysis are data-adaptive decomposition techniques such as singular spectrum analysis (SSA) [16, 17, 18, 19, 20], empirical mode decomposition (EMD) [21, 22] and nonlinear mode decomposition (NMD) [23, 24]. The basic idea of these data-adaptive techniques is to decompose a time series as a sum of modes that describe separately non-oscillating trend, (quasi-)periodic components and high-frequency noise. These techniques are nonparametric because, in contrast to the classical Fourier decomposition, the modes are not model dependent and do not need to be periodic sine or cosine functions. Instead, the modes are derived from the data itself, and they are not limited to a single time scale or a limited range of scales, but describe the data at all scales present. Recently, we applied SSA to quantify irregular rhythms in actigraphy data in the case of persons that suffer from acute insomnia [25], and SSA was applied as well to study irregular patterns in neural and locomotor activity in hamsters [26], but apart from these studies, data-adaptive time-series methods have not been applied in chronobiology. The lack of accessible specialized software to carry out data-adaptive time-series analysis may be one of the reasons that these techniques, to date, have not been applied to circadian rhythm research; fortunately, several open-source implementations have recently become available in multiple platforms such as Mathematica, MATLAB, R, Python, and so on, for SSA [27, 28, 29, 30], EMD [31, 32, 33, 34] and NMD [35].

Another disadvantage of the cosinor method is that it is unable to measure rhythm fragmentation [2, 36, 37]. In actigraphy, rhythm fragmentation was originally defined as the deterioration of the regular circadian rhythm by the occurrence of daytime naps and/or nocturnal activity episodes. On the other hand, spontaneous moment-to-moment fluctuations are a characteristic property of actigraphy time series in particular and of physiological variables in general, and a moderate level of rhythm fragmentation may be indicative of a healthy physiological capacity to respond to random and unforeseen events at multiple time scales. Spectral analysis and other time-series decomposition techniques are ideal tools to study such rhythm fragmentation in actigraphy and other physiological data because they quantify the relative contribution of different time-series components at different time scales to the total variance of the experimental data [25]. A $1/f$ fractal power law may be an indication of such an optimal level of rhythm fragmentation, because it has been observed empirically over a wide range of ultradian time scales in heart rate time series and actigraphy, whereas adverse conditions such as aging [38], cardiovascular disease [39], dementia [40, 41, 42] and insomnia [25, 43] have been found to correlate with deviations of this power law.

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. [43] (see Figure 1). These actigraphy time series show the number of movements per minute for a total duration of 20,160 minutes ( $7\times 24$ 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.

Figure 1.

Two-week continuous actigraphy time series of number of movements per minute for (a) young female adult A (23 years old) with regular circadian cycle, (b) young male adult B (22 years old) with irregular circadian cycle (c) older male adult C (82 years old) with regular circadian cycle. All subjects are asymptomatic controls. Vertical gridlines indicate midnight. Also shown are estimations using cosinor analysis of the circadian 24-h cycle (continuous curve) around a constant mesor (broken line) (see Eq. (2)). All time series have a length of $N=20160$ data points (corresponding to $7\times 24$ h) (data from the public database of Ref. [43]).

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 $x\left(n\right)$ of $N$ discrete data points:

$$x\left(n\right)=\left\{x_1,x_2,\dots x_N\right\}.$$ E1

The procedure consists of fitting a continuous periodic function $y\left(t\right)$ to time series $x\left(n\right)$ :

$$y\left(t\right)=M+A\text{cos}\left(2\mathit{\pi t}/T+\varphi \right)$$ E2

where $M$ is the average value or mesor around which the function oscillates, $T$ is the period, $A$ is the amplitude and $\varphi$ is the phase which defines the value at which the function begins at the start of the monitoring at $t=0$ . By taking into account trigonometric rules, the cosinor function of Eq. (2) can be rewritten as $y\left(t\right)=M+B\text{cos}\left(2\mathit{\pi t}/T\right)+C\text{sin}\left(2\mathit{\pi t}/T\right)$ , where $B$ and $C$ are amplitudes and where the phase is given implicitly by the superposition of the sine and cosine functions. Function $y\left(t\right)$ is fitted to the data by minimizing the summed square residual errors $e_n^2={\left(x_n-y_n\right)}^2$ for all data points $n=1,2,\dots ,N$ , and the value of period $T$ for which the amplitude $A$ 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 $T\approx 24$ h = 1440 min. Once chosen period $T$ , the other circadian parameters $M,A$ and $\varphi$ are determined as well [8]. Phase $\varphi$ 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 ${\varphi }_0$ , 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 ( ${}^{\circ }$ , taking into account the relation ${360}^{\circ }=24$ h), relative to this reference time. The fitted function may be generalized to include more than one period, $y\left(t\right)=M+{\sum }_k{A}_k\text{cos}\left(2\mathit{\pi t}/T_k+{\varphi }_k\right)$ , where the sum usually runs over a small number $k$ of different periods. In the present case, Figure 2 shows many ultradian ( $T<1440$ min) and infradian periods ( $T>1440$ min) that have nonzero amplitudes, but in the present case, there are no other periods than the circadian period $T\approx 24$ h that are clearly distinguishable from the neighboring values to warrant their inclusion in the model function $y\left(t\right)$ . An important result in cosinor analysis is the coefficient of determination $R^2$ , which compares the variance of the residual errors $e_n$ around the fitted model $y$ to the variance of the time series $x\left(n\right)$ around its average value:

$$R^2=1-\frac{\mathrm{Var}\left(\mathrm{e}\right)}{\mathrm{Var}\left(\mathrm{x}\right)}$$ E3

$$=1-\frac{\sum _{n=1}^N{\left(x_n-y_n\right)}^2}{\sum _{n=1}^N{\left(x_n-\text{mean}\left(x\right)\right)}^2},$$ E4

such that $R^2$ is a measure for the fraction of the variance of the time series that can be explained by model $y\left(t\right)$ .

Figure 2.

Regression analysis of Eq. (2) to the time series of subjects A, B and C of Figure 1. Period $T$ is varied from 1 to 2000 min, and amplitude $A$ (shaded curve), mesor $M$ (black line) and phase $\varphi$ (gray curve) are plotted as a function of $T$ . Local maxima near $T=24$ h = 1440 min are indicated with a vertical line and are located at $T=1430$ min (subject A), $T=1438$ min (subject B) and $T=1440$ min (subject C).

Results for the circadian parameters of the cosinor model function y(t) of Eq. (2) for subjects A, B and C are presented in Figure 1. In the present of one single period $T$ , function $y\left(t\right)$ is periodic, and mesor $M$ and amplitude $A$ are constant and capture the average properties of the time series without the possibility to describe day-to-day variability. If $T\equiv 1440$ min, such as for subject C, also acrophase ${\varphi }_0$ is constant, and there is no variability for ${\varphi }_0$ . If, on the other hand, $T\ne 1440$ min, there is a day-to-day phase advance ( $T<1440$ min) or phase delay ( $T>1440$ min), respectively, as is the case for subjects A and B, and one can calculate an average value and variability measures for ${\varphi }_0$ .

		Cosinor			Fourier filter			DWT			SSA
		A	B	C	A	B	C	A	B	C	A	B	C
Mean (M)	(1/min)	222.39	182.81	104.77	222.3	181.9	104.7	222.3	181.9	104.7	227.8	181.8	106.2
Mean (T)	(min)	1430	1438	1440	1409.7	1444.3	1428.1	1409.3	1446.7	1433.3	1421.4	1396.1	1431.7
Mean (A)	(1/min)	154.88	113.66	109.14	174.7	145.4	112.9	151.4	142.0	108.3	165.2	137.5	112.1
Mean ( ${\varphi }_0$ )	( ${}^{\circ }$ )	230.75	240.75	198.75	229.3	217.0	198.6	218.9	210.9	200.1	226.3	220.6	199.9
Mean ( ${\varphi }_0$ )	(hh:mm)	15:23	16:03	13:15	15:17	14:28	13:14	14:36	14:04	13:20	15:05	14:42	13:20
SD (M)	(1/min)	—	—	—	—	—	—	0.8	1.3	0.4	39.3	49.0	15.4
SD (T)	(min)	—	—	—	116.4	261.9	78.0	127.9	326.5	72.5	105.8	190.0	60.2
SD (A)	(1/min)	—	—	—	108.0	87.7	35.5	35.9	46.9	14.8	53.9	59.8	17.9
SD ( ${\varphi }_0$ )	( ${}^{\circ }$ )	10.46	2.09	0	30.9	72.5	15.1	36.6	68.8	13.7	27.8	73.2	11.1
SD ( ${\varphi }_0$ )	(hh:mm)	00:42	00:08	00:00	02:04	04:50	01:00	2:26	04:35	00:55	01:51	04:53	00:44
$R^2$		0.112	0.058	0.207	0.193	0.141	0.238	0.174	0.129	0.230	0.202	0.161	0.241

Table 1.

Circadian parameters mesor $M$ , circadian period $T$ , amplitude $A$ and acrophase ${\varphi }_0$ for subjects A, B and C according to cosinor analysis, Fourier filter, discrete wavelet transform (DWT) using the Daubechies [4] mother wavelet and singular spectrum analysis (SSA) with parameter $L=1440$ . Presented are average values (mean), standard deviation (SD) and the coefficient of determination $R^2$ .

3. Fourier spectral analysis

Fourier spectral analysis makes the supposition that the fluctuations of time series $x\left(n\right)$ with $n=1,2,\dots ,N$ may be interpreted as the superposition of ${\kappa }_{\max }=N/2$ independent harmonic oscillators, where $2\pi /T_{{\kappa }_{\max }}$ is the Nyquist frequency and where each harmonic oscillator corresponds to a periodic function (see Ref. [44]):

$$y\left(t\right)=\sum _{\kappa =1}^{{\kappa }_{\max }}{A}_{\kappa }\text{cos}\left(2\mathit{\pi t}/T_{\kappa }+{\varphi }_{\kappa }\right)$$ E5

$$=\sum _{\kappa =1}^{{\kappa }_{\max }}{B}_{\kappa }\text{cos}\left(2\mathit{\pi t}/T_{\kappa }\right)+C_{\kappa }\text{sin}\left(2\mathit{\pi t}/T_{\kappa }\right),$$ E6

which establishes a link between the time domain $t$ and the frequency domain $f_{\kappa }=2\mathit{\pi t}/T_{\kappa }$ . Here, the first term $\kappa =1$ corresponds to $f=0$ or $T=\infty$ and is the direct current (DC) term around which the other terms $\kappa >1$ oscillate and may be interpreted as the equivalent of the mesor $M$ of cosinor analysis. One of the most important results of Fourier spectral analysis is the power spectrum, which gives the power $P\left(f_{\kappa }\right)$ as a function of frequency $f_{\kappa }$ . The variance of time series $x$ is given by

$$\mathrm{Var}\left(x\right)=\frac{1}{N}\sum _{n=1}^N{\left(x_n-\text{mean}\left(x\right)\right)}^2,$$ E7

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.

$$\mathrm{Var}\left(x\right)=\frac{1}{\left(N-1\right)}\sum _{\kappa =2}^{{\kappa }_{\max }}P\left(f_{\kappa }\right),$$ E8

which allows us to interpret the power $P\left(f_{\kappa }\right)$ of the component with frequency $f_{\kappa }$ as a partial variance, and we can focus our attention to the components that concentrate most of the variance of the time series $x$ .

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 $P\left(f\right)$ as a function of frequency in linear scale show for all subjects a dominant peak at the circadian frequency of $f\approx 14$ 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 $1/f^{\beta }$ scaling with spectral exponent $\beta \approx 1$ and which for subject C appears to be truncated at about $f={10}^2$ oscillations, below which the power spectrum flattens out. A drawback of the power spectrum representation $P\left(f\right)$ as a function of $f$ 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 $\beta$ . The power spectrum $P\left(f\right)$ can be re-ordered in the shape of a so-called scree diagram or Zipf-type plot $P\left(k\right)$ , where the partial variances have been ordered from the most to the least dominant. The advantage of scree diagrams $P\left(k\right)$ is that much less dispersion is present, such that power laws $1/k^{\gamma }$ can be more easily determined, and scaling properties are conserved between both representations, i.e. $\beta =\gamma$ . The disadvantage is that the exact frequency ordering is lost: low-frequency oscillations tend to cluster to the high-magnitude side of the scree diagram and high-frequency oscillations to the low-magnitude noisy tail, but not all components follow this tendency. Scree diagrams for subjects A, B and C are compared in Figure 4. In absolute values (panel a), partial variances of time-series components are comparable for subjects A and B at all scales, whereas those for subject C are smaller. In relative values (panel b), on the one hand, subject C has a larger and subject B a smaller contribution of the circadian cycle in comparison with subject A, but on the other hand, subject C has less and subject B more rhythm fragmentation over a wide range of ultradian frequency scales. Overall, ultradian components appear to scale according to a $1/f$ power law ( $\gamma \approx 1$ ) for subject A, whereas this power law appears to be slightly broken below $k={10}^2$ for subjects B and C. Fourier spectral analysis (Eq. (5)) allows to filter the spectral components of time series. Coefficients of unwanted spectral contributions can be put zero, $A_{\kappa }=0,$ or $B_{\kappa }=$ $C_{\kappa }=0$ , for selected $\kappa$ , to construct a low-pass, high-pass, bandpass or band-stop filter. Figure 5 shows results for a bandpass filter applied to the time series of subjects A, B and C where arbitrary limiting frequencies $f_{\min }=1$ and $f_{\max }=20$ were chosen (to include frequencies above a period of 12 h to obtain a circadian cycle with a single maximum per day). The resulting circadian fit describes variability in time and in amplitude, and corresponding circadian parameters are listed in Table 1.

Figure 3.

Fourier spectral analysis for subjects A, B and C. Power spectrum $P\left(f\right)$ as a function of frequency $f$ in linear scale where the DC term has been omitted for reasons of visibility (left-hand column), logarithmic scale (middle column) and as a scree diagram $P\left(k\right)$ ordered according to magnitude in logarithmic scale (right-hand column). Frequency units are the number of oscillations during the whole duration of the time series.

Figure 4.

Comparison of the Fourier scree diagrams $P\left(k\right)$ for subjects A (green), B (blue) and C (orange), using (a) absolute values and (b) relative values where the total variance of all time series has been normalized Var( $x$ ) = 1.

Figure 5.

Filtered time series of subjects A, B and C of Figure 1 using a Fourier filter between frequencies $f_{\min }=1$ and $f_{\max }=20$ in units of the number of oscillations for the whole duration of the time series of 2 weeks.

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 [45]. 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 $\mathrm{\Psi }\left(t\right)$ , called “mother wavelet”, and the signal $x\left(n\right)$ itself. The resulting CWT depends on $\mathrm{\Psi }$ (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 $t$ and scaled by a factor $s$ :

$$W\left(t,s\right)={\int }_{-\infty }^{\infty }{\mathrm{\Psi }}^{\ast }\left(\frac{u-t}{s}\right)x\left(u\right)\mathit{du},$$ E9

where the asterisk denotes the complex conjugate, and the result is the scalogram $W\left(t,s\right)$ , which represents at each time $t$ the instantaneous period $T$ (represented by scale $s$ ) and the instantaneous intensity or power $W$ of the signal [14]. 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.

Figure 6.

Mother wavelet basis functions, for use with CWT, are (a) Morlet, (b) DGaussian, (c) Gabor, (d) Mexican hat and (e) Paul and with DWT are (f) Haar, (g) Daubechies [4] and (h) BiorthogonalSpline [3, 5]. Basis functions can be real (continuous curve), or complex, in which case the real part (continuous curve) and the imaginary part (dashed curve) are plotted separately.

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 $T\approx 1440$ min, but the intensity and the instantaneous value of $T$ vary in time, reflecting the time and amplitude variability of the circadian rhythm. At ultradian frequencies, $T<1440$ 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 $W\left(t,s\right)$ 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 $T\approx 1440$ min and the $1/f$ 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.

Figure 7.

CWT analysis of time series of subjects A, B and C of Figure 1 using a Morlet mother wavelet. Scalograms are presented in (a)-(c), where vertical gridlines indicate midnight. Power spectra (time-averaged projections on the period axis) are presented in (d) in absolute values and in (e) in relative values (unit variance).

Figure 8.

CWT analysis of time series of subject A of Figure 1 using different mother wavelets: (a) Gaussian, (b) Gabor, (c) Mexican hat and (d) Paul. Vertical gridlines at midnight.

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. [15] 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, ${\text{log}}_2\left(N\right)={\text{log}}_2\left(N\right)\approx 14$ . 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 $R^2$ for circadian fits for different choices of the DWT mother wavelet. Of all choices tried, Daubechies (4) appears to maximize the coefficient of determination $R^2$ .

Figure 9.

DWT time-series components for subject A using Daubechies [4] mother wavelets. Log ${}_2$ (dim) $\approx 14$ scales have been used in the decomposition, and components are ordered from smallest (scale 1) to largest (scale 14). Components are shown with a common time axis (left-hand column). The 10th scale oscillates approximately once every 24 h and can be identified as the circadian cycle, whereas the smaller scales (1–9) correspond to ultradian rhythms and the larger scales (11–14) to infradian rhythms. The last component, without order number, is the residual after all wavelet components (scales 1–14) have been subtracted and can be interpreted as the mesor.

Figure 10.

Estimation of the circadian cycle for subjects A, B and C according to the DWT analysis with Daubechies [4] mother wavelets, using the 10th of 14 scales. Vertical gridlines at midnight.

Figure 11.

Partial variance (%) carried by the scales 1–14 of the DWT analysis for subjects A, B and C, using Daubechies [4] mother wavelets. The last (largest) partial variance corresponds to the residual after all wavelet components have been subtracted of the time series.

DWT	A	B	C
Haar	0.154	0.101	0.205
Daubechies (1)	0.154	0.101	0.205
Daubechies (2)	0.179	0.132	0.220
Daubechies (4)	0.174	0.129	0.230
Daubechies (6)	0.143	0.104	0.223
BiorthogonalSpline (1,3)	0.158	0.103	0.212
BiorthogonalSpline (2,2)	0.095	0.046	0.185
BiorthogonalSpline (2,6)	0.106	0.057	0.194
BiorthogonalSpline (3,5)	0.150	0.096	0.201
BattleLemarie (2)	0.119	0.045	0.195
BattleLemarie (3)	0.127	0.075	0.212
BattleLemarie (4)	0.152	0.103	0.222
BattleLemarie (15)	0.128	0.079	0.210

Table 2.

Coefficient of determination $R^2$ according to DWT analysis using different types of mother wavelet. As suggested in Ref. [15], the Daubechies (4) basis is the most adequate to describe circadian cycles in actigraphy data, and the values given here are the same as listed in Table 1.

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. [19], whereas a larger and very complete review article is in Ref. [20]. We have discussed the SSA method previously in Ref. [25]. 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 data-adaptive or model-independent because the basis functions are generated from the data itself. It can be shown that the sum of all time-series components is identical to the original time series. Below, a summary is given of the most important outcomes of SSA analysis. When applying SSA to a discrete time series $x\left(n\right)$ with length $n=1,\dots ,N$ (see Eq. (1)), a particular window length $L$ must be chosen as an initial parameter, with $2\le L\le N/2$ , which allows to fix the number of components $r$ into which the time series will be decomposed:

$$x\left(n\right)=\sum _{k=1}^r{\sigma }_k{g}_k\left(n\right),$$ E10

where $g_k\left(n\right)$ is the time-series component, ${\sigma }_k$ is the singular value that serves as weights for the components and $r\le \text{min}\left(K,L\right)$ with $K=N-L+1$ . Only (quasi-)periodicities with average length $T\lesssim L$ will be resolved into separate time-series components, whereas those with lengths $T>L$ will be absorbed in the trend component. One can chose $L$ as a multiple of the (average) periodicity of the data, i.e. $L=\mathit{mT}$ , where $m$ is an integer number. In the case of circadian data, the obvious choice would be $L=T=24$ h = 1140 min. It can be shown that in the limit for $L\to N/2$ , SSA converges to Fourier spectral analysis [20], where a time series is always decomposed as the superposition of $N/2$ -independent oscillators (Nyquist theorem). Whereas the Fourier power spectrum of a quasiperiodic time series with average period $T$ would correspond to a broad Gaussian peak around the central frequency $f=1/T$ , intuitively, it can be understood that for an adequate choice of the parameter $L$ the neighboring Fourier components of this broad Gaussian peak can be “compressed” within a single SSA component. One of the main results of SSA analysis is the so-called scree diagram that visually represents the partial variances ${\lambda }_k={\sigma }_k^2$ , ordered according to magnitude from the most to the least dominant, where ${\lambda }_k$ can be interpreted as the variance of the “subphase space” of time-series component $g_k\left(n\right)$ and where ${\lambda }_{\mathit{tot}}={\sum }_{k=1}^r{\lambda }_k$ is the total variance of the phase space of the original time series $x\left(n\right)$ . The dominant partial variance ${\lambda }_1$ , associated with component $g_1\left(n\right)$ , usually corresponds to the trend. Dominant periodicities can be recognized as “steps” in the scree diagram, i.e. two successive partial variances ${\lambda }_k$ and ${\lambda }_{k+1}$ that are nearly degenerated and clearly distinguishable from the neighboring partial variances, where the corresponding components $g_k\left(n\right)$ and $g_{k+1}\left(n\right)$ are the Fourier equivalents of a sine and cosine function with the same frequency. Higher-order partial variances ${\lambda }_k$ tend to have values that decrease gradually and continuously with $k$ , indicating that at these scales it is impossible to distinguish any individual time-series components $g_k\left(n\right)$ that separately can be assigned physical significance.

Figure 12 shows some details of the decomposition of the time series of subject A. In the scree diagram, a trend mode $g_1\left(t\right)$ , a mode $g_2\left(t\right)+g_3\left(t\right)$ which we will identify as the circadian (24 h) cycle and a mode $g_4\left(t\right)+g_5\left(t\right)$ 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 $1/k^{\gamma }$ with $\gamma \approx 1$ . The weighted correlation matrix shows that the trend mode and the circadian mode are uncorrelated from the higher-order modes, but the $\approx 12$ -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 $\approx 12$ -h ultradian mode. Figure 13 shows the average correlation of the circadian mode with other SSA components, and broad minima can be observed for $L$ being equal or a multiple of the circadian periodicity, $L=m\times 1440$ 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 $L=1140$ 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 ( ${10}^{1.5}\le k\le {10}^{3.0}$ ), all subjects show a very similar behavior with a $1/f$ scaling. At the lowest ultradian frequencies ( ${10}^{0.5}\le k\le {10}^{1.5}$ ), $1/f$ scaling continues for subject A, whereas rhythm fragmentation is increased for subject B and decreased for subject C.

Figure 12.

SSA analysis of time series of subject A of Figure 1. Shown are (a) the scree diagram of partial variances ${\lambda }_k$ , (b) the matrix of weighted correlations between the first 20 time-series components $g_k\left(t\right)$ and $g_l\left(t\right)$ with $k,l=1,\dots ,20$ and (c) the first five time-series components $g_k\left(t\right)$ with $k=1,\dots ,5$ , of which $g_1\left(t\right)$ is the trend component, $g_2\left(t\right)+g_3\left(t\right)$ is the circadian rhythm and $g_4\left(t\right)+g_5\left(t\right)$ is the 12 h ultradian rhythm.

Figure 13.

Average correlation the circadian mode $g_2\left(t\right)+g_3\left(t\right)$ with the neighboring components $g_1\left(t\right),g_4\left(t\right),\dots g_{20}\left(t\right)$ as a function of parameter $L$ for subjects. Average correlation is lower for subject C (orange curve) than for subject A (green curve), which is lower than for subject B (blue curve). Gridlines indicate multiples of the 24-h period, 1440 and 2880 min.

Figure 14.

Estimation of the circadian cycle of subjects A-C according to SSA analysis using $L=1440$ .

Figure 15.

Scree diagram of SSA analysis of subjects A-C, showing (a) partial variances and (b) fractional partial variances.

6. Discussion

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 R² of Table 1, the Fourier filter, DWT and SSA describe the circadian cycle better than the cosinor analysis with one single period, and SSA gives the best description of all methods discussed here. One of the reasons may be that cosinor cannot take into account the variability in time and amplitude of the experimental time series, whereas the other methods can. It is less clear why SSA analysis results in the best fit, the average amplitude and the variability of the circadian parameters tend to be larger for the Fourier filter than for SSA, but they tend to be smaller for DWT. The goodness of fit to the data for DWT depends on the specific mother wavelet used, but we chose the Daubechies (4) mother wavelet because of the maximized $R^2$ for all the different mother wavelets that we experimented with. On the other hand, the number of DWT scales might be increased or decreased, in order to adjust the variability of the circadian mode for a better fit to the experimental data, but there is no rule of thumb that says how many scales to choose. It is possible that the Fourier filter description may be improved by carefully adjusting $f_{\min }$ and $f_{\max }$ , but it is not a priori clear which limiting frequencies will result in the best $R^2$ . Figure 13 suggests that it might be possible to slightly improve the SSA description by fine tuning the parameter value $L$ to a value for which a global minimum is obtained in the correlation of the circadian mode $g_2\left(t\right)+g_3\left(t\right)$ , but on the other hand, the broad minima basins suggest that the calculations are rather stable. Thus, as long as $L\approx 1440$ , the precise value of $L$ and the number of components into which the time series is decomposed have little influence on the description of the circadian cycle by SSA. Table 3 compares the similarity of the description of the circadian cycle between the different methodologies. Results for the Pearson product-moment correlation $r$ and Spearman’s rank correlation $\rho$ are very similar. When the circadian cycle is very regular, as for subject C, descriptions by different methods resemble very much. On the other hand, when the circadian cycle is very irregular, as for subject B, the different methods of cosinor, Fourier filter, DWT and SSA give different results.

A		r				B		r				C		r
		cosinor	filter	DWT	SSA			cosinor	filter	DWT	SSA			cosinor	filter	DWT	SSA
$\mathit{\rho }$	cosinor		0.758	0.845	0.840	$\mathit{\rho }$	cosinor		0.646	0.757	0.791	$\mathit{\rho }$	cosinor		0.933	0.968	0.976
	filter	0.793		0.881	0.949		filter	0.638		0.830	0.906		filter	0.930		0.964	0.977
	DWT	0.843	0.912		0.950		DWT	0.758	0.820		0.949		DWT	0.960	0.963		0.992
	SSA	0.842	0.953	0.950			SSA	0.790	0.894	0.947			SSA	0.968	0.976	0.991

Table 3.

Correlation between the calculation of the circadian cycle according to the cosinor method, the Fourier filter, DWT and SSA, for subjects A, B and C. Results are given for the Pearson product-moment correlation $r$ (upper-right triangle) and Spearman’s rank correlation $\rho$ (bottom-left triangle).

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 $1/f$ 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.

7. Conclusions

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.

Acknowledgments

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.