Circadian parameters mesor

## 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.

### Keywords

- singular spectrum analysis
- SSA
- wavelets
- spectral analysis
- Fourier analysis
- data-adaptive
- model-free

## 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

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 (

## 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

The procedure consists of fitting a continuous periodic function

where

such that

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

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 ( | ( | 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 ( | (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 ( | ( | 10.46 | 2.09 | 0 | 30.9 | 72.5 | 15.1 | 36.6 | 68.8 | 13.7 | 27.8 | 73.2 | 11.1 |

SD ( | (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 |

0.112 | 0.058 | 0.207 | 0.193 | 0.141 | 0.238 | 0.174 | 0.129 | 0.230 | 0.202 | 0.161 | 0.241 |

## 3. Fourier spectral analysis

Fourier spectral analysis makes the supposition that the fluctuations of time series

which establishes a link between the time domain

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 * partial variance*, and we can focus our attention to the components that concentrate most of the variance of the time series

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 * scree diagram*or

*-type plot*Zipf

## 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

where the asterisk denotes the complex conjugate, and the result is the scalogram

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

### 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,

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 |

## 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

*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*model-independent

where * singular value*that serves as weights for the components and

partial variances

Figure 12 shows some details of the decomposition of the time series of subject A. In the scree diagram, a trend mode

## 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 ^{2} 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

A | B | C | |||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

cosinor | filter | DWT | SSA | cosinor | filter | DWT | SSA | cosinor | filter | DWT | SSA | ||||||

0.758 | 0.845 | 0.840 | 0.646 | 0.757 | 0.791 | 0.933 | 0.968 | 0.976 | |||||||||

0.793 | 0.881 | 0.949 | 0.638 | 0.830 | 0.906 | 0.930 | 0.964 | 0.977 | |||||||||

0.843 | 0.912 | 0.950 | 0.758 | 0.820 | 0.949 | 0.960 | 0.963 | 0.992 | |||||||||

0.842 | 0.953 | 0.950 | 0.790 | 0.894 | 0.947 | 0.968 | 0.976 | 0.991 |

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

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.