Circadian rhythms are ∼ 24 h oscillatory physiological processes [1]. The rhythms are synchronized with external cues, such as day-night cycles, by the endogenous central pacemaker, called the circadian clock, which is in the hypothalamic suprachiasmatic nucleus (SCN) [1]. Disruption of this synchrony, for instance, due to a shift-work lifestyle, increases the risk for many diseases, such as psychiatric disorders, neurodegenerative diseases and cancer [2,3]. Thus, identifying optimal sleep schedules for shift workers that lead to minimizing fatigue has received attention [4]. Moreover, optimal drug and surgery timing has been considered an inevitable part of precision medicine [5]. Notably, timed infusions of anti-cancer drugs can improve the survival of cancer patients [6,7]. These clinical applications based on circadian physiology require accurate identification of underlying circadian rhythms in physiological processes [8].
Previously, several methods to track the clock state have been developed [9,10]. For example, one of the widely used approaches is to measure the secretion onset time of the sleep-regulating hormone, melatonin, from saliva or blood samples under dim-light conditions [9]. Although the assessment of dim-light melatonin onset (DLMO) is well validated and thus considered the gold standard, challenges arise when using DLMO to develop chronotherapeutics [5,8]. Major complications include: (i) DLMO assessment is labour-intensive and invasive [11]; (ii) the sample collection occurs in laboratory settings, making it infeasible to measure DLMO continuously for a long time (e.g. greater than 30 days) [11]; (iii) the subjective choice of threshold for DLMO assessment (e.g. 3 or 4 pg ml -1 ) can affect the clock state estimation [12]; and (iv) it is costly [13]. Thus, large-scale epidemiological studies with DLMO are challenging.
One possible remedy for this is to exploit wearable technologies. Wearables can continuously measure behavioural and physiological signals such as activity, heart rate (HR) and body temperature (BT) over long periods in real-world settings (figure 1a) [8,16]. These signals are influenced by the circadian clock and thus have circadian variation [17,18]. Accordingly, tracking of the underlying circadian rhythms from wearable measurements enables real-world chronomedicine. However, wearable technologies pose their own set of challenges. The HR and BT circadian rhythms can become obscured by many factors, such as hormones, sleep, activity, caffeine and stress [19]. The obscured circadian signals are further hidden in the measurement noise [14]. As a result, only noisy daily (i.e. time-of-day) rhythms that arise from the combined influence of circadian and state-evoked external components are typically observed.
To overcome the challenges and extract circadian information specifically related to endogenous timers, several mathematical and statistical methodologies have been proposed [14,16,20]. Milestone studies done by Brown and Czeisler showed that the endogenously generated BT rhythms with a period 24 h (i.e. the circadian rhythm in BT) can be identified from constant-routine core body temperature data collected in laboratory settings using harmonic-regression models [15]. Other important early work showed that the characteristics of the circadian rhythms in BT, such as the type-0 phase response curve to a three-cycle bright light stimulus and the timing of core body temperature minimum, can be accurately captured by van der Pol type oscillator models [21]. Based on these foundational studies, researchers have recently made meaningful progress in identifying circadian rhythms from wearable data. The dashed grey lines represent the true values of the parameters described in electronic supplementary material, table S1. (e,f ) Box-and-whisker plots of 10 2 mean estimates that were obtained when increasing the measurement interval. The estimates of BT (e) and HR ( f ) parameters are still reasonably accurate even if the measurement interval is increased. Here, the two-harmonic model and the single harmonic model were used for the estimation of the BT and HR parameters, respectively, as done in the previous studies [14,15].
Huang and colleagues showed that a mathematical model taking wearable activity data can predict the DLMO in normal free-living conditions [20]. Moreover, a statistical framework has been recently proposed to extract the circadian rhythm in HR that originates in the sinoatrial node of the heart [14]. This method has been tested against HR data from a constant-routine protocol, and HR phase estimates from this method display different dynamics than sleep or activityderived timings. The HR phase estimates generally aligned with DLMO predictions during normally entrained scenarios, but became desynchronized during social distancing [22]. Together, these findings suggest an ability to capture endogenous circadian signals from algorithms applied to wearable data. Considering this, we will consistently employ the term 'circadian' in the rest of the paper to refer to the 24 h rhythm signals extracted using filtering frameworks. Despite this progress, there are still apparent challenges. Specifically, it is challenging to use simple and efficient methods such as least-squares methods (LSMs) for the estimation of circadian parameters because the necessary assumption that the noise process is independent Gaussian does not hold in physiological time-series data. Indeed, external and internal factors lead to autoregressive noise processes in BT and HR [14,15]. Due to this, previously proposed methods are based on a more computationally expensive Bayesian inference framework with Markov chain Monte Carlo (MCMC) [14,15]. As a result, these methods cannot perform circadian assessment solely by the computing power of consumer-grade wearables. Thus, time-consuming intermediate steps are required, such as the anonymous transmission of wearable data to secure computer clusters. Moreover, the identification of a suitable model to analyse wearable data remains to be studied. While it has been found that the two-harmonic-regression model is a reasonable choice to extract circadian information from BT data measured in a constant routine protocol [15], the suitability of this model for wearable data needs to be explored. For this, a systematic comparison of the goodness-of-fit measures of candidate harmonic models applied to large-scale wearable datasets is required, which is impractical with previous computationally demanding methods [14,15].
In this study, we developed an efficient method to extract physiological parameters from wearable data. It transforms a harmonic-regression model with correlated noise into that with independent Gaussian noise under the assumption that the measurement interval is sufficiently small. This allows us to use LSMs to efficiently estimate the circadian parameters and their uncertainty. We showed that our method is reasonably accurate and computationally much more efficient than the previous methods [14,15] by testing it on in silico BT and HR data. Using our method, we also found that the time of the circadian minimum, referred to as the circadian phase in previous studies [14,15], can be more accurately estimated with a single harmonic model than with a multiple harmonic model having more degrees of freedom. We next applied our method to real-world BT and HR data: ∼ 600 days of BT data and 184 000 days of HR data, revealing the inter-and intra-individual differences in real-world dynamics of circadian rhythms. Finally, we combined this method with a Kalman filter (KF) framework [23], which enables efficient real-world tracking of the state of the molecular clocks in tissues. The computational gains, systematic model testing, and real-world applicability of this work set the stage for personalized digital medicine with wearable devices.
2. Methods 2.1. Body temperature parameter estimation from wearable data 2.1.1. A harmonic-regression model of the human body temperature rhythm
The human BT oscillates with ∼ 24 h periodicity [24]. This circadian oscillation has been analysed with Fourier series representations (i.e. harmonic-regression models) [15,16,[25][26][27]. Specifically, Brown and Czeisler identified that the first-order autoregressive noise process needs to be incorporated into the harmonic models to successfully analyse the BT [15]. Following this previous study, we adopted a harmonic-regression-plusfirst-order-autoregressive model (equation (2.1)) to analyse our wearable BT data collected from hospital cancer patients who are unable to perform an excessive physiological activity that may affect the BT (electronic supplementary material).
where
and it is assumed that the circadian period τ = 24 h, |α| < 1, and the e t values are distributed as Gaussian random variables with mean zero and variance s 2 e . The autoregressive noise process describes the ongoing effects of external factors on the BT. The Gaussian noise e t represents new external influences, for example, from stress, hormones and measurement error. The autocorrelation factor α represents the ongoing contribution of external factors.
To efficiently estimate the parameters of the BT model (equation (2.1)), we propose an approximation-based linear LSM. We first reformulate and approximate equation (2.1) as follows:
e t : ð2:2Þ The second equality in equation (2.2) holds as v t is a first-order autoregressive noise process (equation (2.1)). In the fourth line in equation (2.2), the assumption (s t ≈ s t-Δt ) is used, which is reasonable if the measurement interval Δt is small. Typically, Δt is sufficiently small because wearable devices continuously collect physiological signals with high-frequency resolution (e.g. less than 6 min interval) [8,14,16]. Because the approximated equation (equation (2.2)) only includes independent Gaussian noise, we can exploit the standard linear LSM [28,29] with the approximated equation to estimate the parameters from BT data. Specifically, let y = (y Δt , y 2Δt , …, y t-Δt , y t ) 0 , u ¼ ðm, ã1 , b1 , . . . , ãn , bn , aÞ 0 and e ¼ ðe Dt , e 2Dt , . . . , e tÀDt , e t Þ 0 where m ¼ ð1 À aÞ Á m, ãr ¼ ð1 À aÞ Á a r , and br ¼ ð1 À aÞ Á b r . Then, equation (2.2) can be rewritten in matrix notation as y ¼ Xu þ e, ð2:3Þ where X ¼ 1 cos 2pDt t À Á sin 2pDt t À Á Á Á Á cos 2pnDt t À Á sin 2pnDt t À Á y 0 . . . . . . . . . . . . . . . . . . . . . 1 cos 2pt t À Á sin 2pt t À Á Á Á Á cos 2pnt t À Á sin 2pnt t À Á y tÀDt 2 6 6 4 3 7 7 5 : By applying the linear LSM to equation (2.3), we can derive the probability density function for û that minimizes | |y -Xθ| | given y and X, where | | • | | denotes the 2-norm (equation (2.4)).
û NððX 0 Á XÞ À1 Á X 0 Á y, ðX 0 Á XÞ À1 Á h 2 Þ, ð2:4Þ
where
and dimðÁÞ represents the dimension of a vector. Then, we can compute the probability density of parameters μ, a r , b r and α from the distribution of û (equation (2.4)) using the equations m r ¼ mr =ð1 À aÞ, a r ¼ ãr =ð1 À aÞ and b r ¼ br =ð1 À aÞ. From these, we can compute the mean and uncertainty estimates of the phase and amplitude of the signal s t (electronic supplementary material).
The human HR shows a circadian variation [17]. This HR circadian signal has also been successfully analysed with a harmonicregression-plus-first-order-autoregressive model [14,22,30], as in the analysis of the BT rhythm. One difference between the HR model and the BT model is that the HR model contains a term describing the increase in HR by activity, specified as follows:
where a t is the activity level (i.e. step count) at time t and d is the increase in HR per step (HRpS).
We next describe the nonlinear version of the approximationbased LSM (ALSM) to efficiently estimate the HR model parameters. Like the linear version, we first reformulate and approximate equation (2.5) as follows:
ð2:6Þ
Due to the reformulation and approximation (equation (2.6)), the noise process is converted from autoregressive noise to independent Gaussian noise. Thus, we can now calculate the mean estimate of u ¼ ðm, ã1 , . . . , ãn , b1 , . . . , bn , d, aÞ 0 , denoted by u, and its uncertainty (i.e. the covariance matrix S u ) by using any LSMs for nonlinear squares fitting problems. In this study, we used the Levenberg-Marquardt algorithm (see electronic supplementary material and [31]).
To check the accuracy and precision of the ALSM, we first tested it on in silico wearable data. Specifically, we generated in silico BT and HR data by simulating the previously proposed BT and HR models [14,15] (see electronic supplementary material, table S1). Then, the generated BT and HR data were analysed by the ALSM, which extracted the underlying circadian rhythms in the data and the effect of activity on HR (figure 1b and electronic supplementary material, figure S1).
As a result, a linear version of ALSM can estimate the probability distributions of four BT parameters: mesor, half the range of the fitted oscillatory signal (i.e. amplitude), time of BT minimum (i.e. phase), and autocorrelation strength describing the effect of confounding factors (e.g. stress) on BT (figure 1c, see Methods). A nonlinear version of ALSM can estimate the probability distributions of five HR parameters: mesor, amplitude, phase, autocorrelation strength and the increase in HRpS (figure 1d, see Methods). Note that the presence of red fuzz at the far right edge of the centre in figure 1d does not indicate that the phase estimate samples significantly differ from the true phase. This is because the phase estimates are cyclic quantities measured in units of time of day. Importantly, the estimation of BT and HR parameters remained reasonably successful even when the measurement interval was extended to 6 min, which exceeds the typical sampling interval of wearables (figure 1e,f ). Even with further extension of the measurement interval, the successful estimation of BT parameters persisted (electronic supplementary material, figure S2A). However, the accuracy of estimating the HR phase and amplitude decreased rapidly (electronic supplementary material, figure S2B). To understand the differential effects of increasing the measurement interval on the estimation accuracy of our method for HR data compared with BT data, we examined the relationship between the error resulting from the approximation (s t ≈ s t-Δt ) and both the measurement interval and the model parameters (electronic supplementary material). We found that the upper limit of the error is governed by ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi
where Δt is the measurement interval and a r and b r denote the harmonic coefficients. This suggests that as the amplitude of the rhythmic signal increases, the resulting approximation error also increases. This explains why the impact of increasing the measurement interval on estimation accuracy was more significant in HR data compared with BT data: since the amplitude of the HR rhythm was greater than that of the BT rhythm, the possible approximation error was larger in HR data compared with BT data (electronic supplementary material, figure S2C). Therefore, with an increase in the measurement interval, the estimation accuracy in HR data deteriorated at a faster rate compared with BT data.
Previous experimental studies showed that cardiac rhythmicity is regulated differently during sleep [19]. Thus, the previous algorithm to extract physiological parameters from wearable HR data does not use the data collected during sleep [14,22,30]. Specifically, the model is fitted to 2-day intervals centred at the sleep episode in between. We checked whether the ALSM is still accurate when the data of 2-day intervals except for the data in the sleep period are only given for estimation, as in the previous studies [14,22,30]. Indeed, the HR estimates obtained with the ALSM using the 2-day data collected during the wake period were accurate (figure 2a,b). In particular, even when the measurement interval was increased to 6 min, the estimates of the ALSM were still accurate and precise (figure 2c). This indicates that our method can be used to explore real-world dynamics of the circadian rhythm from the wearable measurements in the same way as the previous method [14].
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 20: 20230030
Our algorithm can identify the underlying circadian rhythms hidden in the influence of confounding factors (e.g. activity) and noise (figures 1 and 2a-c). This provides an opportunity for tracking circadian rhythms in a consecutive manner. We tested this possibility by applying our algorithm to in silico data mimicking typical human lifestyles. First, we simulated a regular real-world lifestyle of humans for 60 days. The HR and activity were simulated as described in electronic supplementary material, except that randomness was introduced into sleep offset and onset time. Specifically, the sleep offset time t off i and onset time t on i of the day i were defined to be t off
t Þ where σ t denotes the randomness of sleep onset and offset times (electronic supplementary material, Table S1). The double-plotted actograms generated by the activity data are shown in figure 2d, with the true timing of the HR minimum set as 3.00. Indeed, the true HR phase was successfully tracked by the ALSM despite the additional randomness in activity data (figure 2d). More importantly, when the activity pattern is suddenly delayed or advanced, resulting in the misalignment between the activity and the HR observed in the previous study [14], the ALSM can successfully track the change of the HR Here, the single harmonic function was used as the HR model as done in the previous study [14]. The box-and-whisker plots were obtained as in figure 1f. (d,e) Double-plotted actograms for virtual subjects generated from in silico activity data. The mean and standard deviation of the HR phase estimates are represented as a red line and range, respectively, with daily activity patterns (black). The ALSM can track the HR phase when the subject is under a constant routine (d ). Moreover, even if the subject has a shift-work lifestyle resulting in the misalignment between the HR and the activity, the HR phase can be estimated (e). ( f ) The ALSM is more computationally efficient than the recently proposed Bayesian MCMC method [14]. The averaged computational time was obtained by applying the two methods to the data 10 times.
phase (figure 2e). This indicates that our method can be used for studying real-world dynamics of the circadian rhythm in HR. The method recently proposed by Bowman and colleagues [14] can also track the HR circadian rhythm. To demonstrate the benefit of the ALSM compared with the previous method, we compared their computational cost. Remarkably, our method was ∼ 300-fold faster than the previous method (figure 2f ). For instance, to analyse data collected for 100 days, the previous method needed ∼ 13 min, while our method needed ∼ 2 s. Such a large improvement in computational efficiency is mainly because our method computes parameter estimates in a deterministic manner while the previous one is based on the Bayesian MCMC framework. Note that we only collected 500 posterior samples using the MCMC when calculating the computational time of the previous method. If more posterior samples need to be collected to improve estimation accuracy, the previous method might become more computationally intractable. Considering this, it is not suitable for large-scale epidemiological studies. In this circumstance, our method can be a promising alternative.
wearable data more accurately with a single-harmonic model than with a multiple-harmonic model
To accurately estimate the circadian phase from wearable data with our method, a suitable harmonic curve-fitting model needs to be adopted. To identify the model, we compared the performance of our method with harmonic-regression models having a different number of harmonics. Specifically, we first generated in silico HR data as in figure 2a, except that the 12 h rhythm was introduced using a two-harmonic model (see electronic supplementary material and table S1 for details). Then, we extracted the time of minimum of the 24 h rhythm (i.e. circadian phase in HR) from the data using our method with either a single harmonic or two harmonics (figure 3a). Interestingly, the phase estimated with the single harmonic matched the true phase more closely than that with the two harmonics, even though the fitted data were generated with the two harmonics. We further compared the performance in the presence of gaps and random data loss, which typically occur in real-world wearable data. Specifically, we generated 10 3 versions of in silico HR data having random circadian phase values drawn from a uniform distribution between 0 and 24 h with multiple-harmonic models (figure 3b,c). We then introduced two varieties of missing data: (i) a gap of 8 h starting at a random point in the data, repeated for each 24 h period and (ii) random data loss with probability 5%. These two sources combined to simulate missing data due to device charging (e.g. during sleep) or collection errors. The generated data were fitted using our method with either a single-harmonic model or a multiple-harmonic model (figure 3b,c). This showed that the single-harmonic model only considering 24 h rhythms outperformed the multiple-harmonic models overall, although the number of harmonics mismatched that in the models used to generate the fitted data.
In figure 3d,e, we quantified the differences between the true phase and the estimate. If the number of harmonics used to generate the data matched that used for estimation, the most accurate estimation of the circadian phase in HR occurred when the number of harmonics (n) was equal to 1 (average error (AE) = 0.671), followed by n = 2 (AE = 1.212) and n = 3 (AE = 2.376) (figure 3d). Even if there was a mismatch between the number of harmonics in the data and the number of harmonics in the model, the single (n = 1) harmonic was still the best model to estimate the circadian phase (figure 3e,f ) as illustrated in figure 3b,c. Specifically, the combination of n = 2 harmonic data and n = 1 harmonic model yielded more accurate predictions than the n = 2 harmonic model and data (AE ¼ 1:012 versus 1:212, figure 3e). Similarly, more accurate phase estimation occurred with the n = 3 harmonic data and n = 1 harmonic model than with the n = 3 harmonic model and data (AE ¼ 1:063 versus 2:376, figure 3f ). Even when the measurement interval was increased and the length of the gap varies, the singleharmonic model overall outperformed the multiple-harmonic models (electronic supplementary material, table S2). Taken together, the single-harmonic model was a suitable choice to estimate the circadian phase in HR even in the presence of multiple harmonic components (i.e. ultradian rhythms) in given wearable data.
Interactions between the amount of data, the measurement interval, and the harmonic order also impacted the estimation of other parameters than circadian phase. Importantly, our method captured the basal HR, increase in HR due to one step, and autocorrelated noise parameters with high accuracy when applied to in silico data generated while varying the gap length between 4, 6 and 8 h, and varying the measurement interval between 1, 2, 4 and 6 min (electronic supplementary material, table S3). The amplitude parameter was captured less accurately overall, and was particularly overestimated for higher numbers of harmonics in the presence of large gaps. This is probably due to the higher-order harmonic terms in the model no longer being independent from the 24 h component due to missing data, and this further supports the performance of the single harmonic mode in capturing parameters in these settings.
We now tested our method on the axillary temperature measurements previously collected using a FDA-approved smart wearable thermometer (TempTraq) [16]. Seventy-two hospital patients with cancer who received chimeric antigen receptor T-cell or haematopoietic cell transplant wore a wearable axillary skin patch to monitor BT every 2 min (see electronic supplementary material). This allowed for the collection of 598 days of densely sampled wearable BT data under non-fever conditions, which provides an opportunity to study the BT dynamics of cancer patients in real-world settings.
We first analysed the high-frequency BT measurements with the ALSM by varying the number of harmonics from one to three (figure 4a). Then, we considered the following four widely used goodness-of-fit statistics: the root-meansquare error (RMSE), the coefficient of determination denoted by R 2 , the Akaike information criterion (AIC), and the Bayesian information criterion (BIC) (see electronic supplementary material, [15,28,32]). Of all the models we considered, the zero-order harmonic that does not consider oscillatory variations in BT described the BT data least well. Specifically, RMSE, AIC and BIC decreased, and R 2 increased when the BT rhythmicity was considered in the parameter estimation (i.e. when the number of harmonics, n > 0) (figure 4b-e and table 1). Taking this into account, oscillatory models are appropriate to describe the key BT wearable data features.
We found that increasing the number of harmonics can improve goodness-of-fit (figure 4a-e and table 1). However, when increasing n = 2 to 3, there were no remarkable improvements in graphical goodness-of-fit (figure 4a) and the statistics (figure 4b-e). This suggests that the twoharmonic model with the autoregressive noise process is a reasonable model for extracting hidden physiological parameters from the wearable BT data. Interestingly, our results were consistent with the model selection findings with BT data, which are typically measured with a rectal thermistor in a constant routine protocol that reduces the effects of external factors and makes endogenous circadian signals more easily observable [15]. Specifically, Brown and Czeisler found that the three-harmonic model with the first-order autoregressive noise process is the best statistical model to describe constant-routine core BT data in terms of the goodness-of-fit measures. However, they also showed no significant improvements in fitting results between the two and three harmonics, resulting in their conclusion that the twoharmonic model is reasonable. Recently, the chest surface body temperature measured by wearables and the core temperature recorded using electronic ingestible pills with radio-frequency transmissions were scrutinized using spectral analysis. This showed that both the chest BT data and core BT data exhibit a dominant period of 24 or 12 h [27], which provides additional support to our results.
We next analysed the fitted BT parameter values from 598 days by varying the number of harmonics (figure 4f-i). The mesor estimates did not change when increasing the number of harmonics (figure 4f and table 1). The mesor estimates fell between 34 C and 37 C, which closely matched the normal range of the temperature of the skin surface of the trunk (33:5 CÀ36:9 C) reported in experimental literature [33,34]. The estimates of amplitude defined to be half the range of the oscillatory signal increased when increasing the number of harmonics (figure 4g and table 1). It is noteworthy that our amplitude estimate obtained with the reasonable two-harmonic model (0:896 C) was consistent with the previously reported one (0:88 C) obtained by fitting the two harmonics to the BT data collected using a wearable chest patch (Movisens, Karlsruhe, Germany) [27]. Moreover, the amplitude estimates of 24 h rhythm (i.e. circadian amplitude) calculated from the estimates of the first-order harmonic coefficients were 0:6 C regardless of the number of harmonics used for estimation (electronic
38.0 35.0 0 1 2 36.5 (b) (a) 24 3 -3 0 RMSE (ºC) R 2 the number of harmonics ( f ) BIC temp. (ºC) 0.8 0 0.4 mesor (ºC) 32 36 40 amp. (ºC) 0 1 . 5 3 . 0 density 2 0 1 phase (h) 0 12 24 0.12 0 0.06 autocorr. strength 0.6 0.8 1.0 20 0 10 0 2 1 3 the number of harmonics 0 1 2 clock time (h) 24 0 12 24 0 12 24 ×10 3 3 -3 0 AIC ×10 3 -1 1 0 0 1 2 3 0 2 1 0 1 2 3 0 1 2 3 0 1 2 3 (c) (g) (d) (h) (e) (i) 0 2 1 3 the number of harmonics and BIC (e) that were computed for 598 days of wearable BT data by varying the number of harmonics used for estimation. ( f-i) Smooth kernel histograms for the probability density function of the mean estimates of mesor ( f ), amplitude (g), phase (h) and autocorrelation strength (i) for 598 days of the BT data. See table 1 for the mean and the standard deviation of the estimates.
Table 1. Goodness-of-fit measures for the harmonic-regression models of BT and the estimated parameters. Here, we analysed the high-frequency BT data collected using a wearable sensor in [16]. The mean of either the measures or the estimates, and their standard deviation (s.d.) are represented in the form of 'mean ± s.d.' in the table. The uncertainty estimate denotes the standard deviation of the estimated probability density of the parameters.
The number of harmonics, n n= 0 n = 1 n = 2 n = 3 goodness of fit measures RMSE 0.811 ± 0.320 0.692 ± 0.261 0.647 ± 0.245 0.641 ± 0.234 R 2 -0.003 ± 0.012 0.241 ± 0.175 0.329 ± 0.182 0.389 ± 0.184 AIC -410.40 ± 590.21 -628.80 ± 590.29 -723.85 ± 590.56 -797.12 ± 595.11 BIC -411.38 ± 590.21 -629.76 ± 590.29 -724.80 ± 590.56 -798.06 ± 595.12 mean estimate mesor ( C) 35.983 ± 0.639 35.985 ± 0.636 35.985 ± 0.636 35.985 ± 0.636 amp. ( C) 0.646 ± 0.373 0.896 ± 0.449 1.068 ± 0.504 phase (h) 10.183 ± 4.179 10.256 ± 4.174 10.162 ± 4.234 autocorr. strength 0.940 ± 0.059 0.921 ± 0.066 0.909 ± 0.071 0.899 ± 0.076 uncertainty estimate mesor ( C) 0.240 ± 0.142 0.168 ± 0.092 0.142 ± 0.078 0.124 ± 0.067 amp. ( C) 0.218 ± 0.118 0.228 ± 0.129 0.230 ± 0.134 phase (h) 2.150 ± 1.406 2.575 ± 1.836 2.498 ± 1.743 autocorr. strength 0.012 ± 0.007 0.014 ± 0.007 0.015 ± 0.008 0.016 ± 0.008 supplementary material, table S4), which closely matched the previously reported value (approx. 0:5 C) obtained using the BT data measured in controlled laboratory settings [18,35]. All possible estimates of phase defined to be the time of skin BT minimum were observed (figure 4h) because we studied a population of cancer patients typically having disrupted BT rhythms [2,36]. The population mean phase in skin BT was ∼ 10.00 even when varying the number of harmonics (figure 4h and table 1), which was consistent with that reported in the previous studies [37,38]. The estimates of autocorrelation strength were greater than or equal to ∼ 0.9, although they slightly decreased when increasing the number of harmonics (figure 4i and table 1). This indicates that a correlated-noise structure needs to be considered to accurately describe non-circadian fluctuations in wearable BT data, resulting in reliable extraction of circadian physiology. The average of the autocorrelation strength estimates with the two harmonics was 0.909, leading to a correlation time of ∼ 1 h. This value was nearly identical to the previously reported estimate (0.897) obtained using constant-routine core BT data [15]. Taken together, most of the circadian and non-circadian parameters extracted from wearable BT data using our method matched the parameters measured in controlled laboratory conditions. Our algorithm outputs both mean estimates and error estimates that describe the statistical uncertainty in the estimates. Specifically, the standard deviation of the probability density associated with the parameters, denoted by the 'uncertainty estimate' below, can be computed (table 1; electronic supplementary material, figure S3 and table S4). The uncertainty estimates typically act as a measure for the reliability of the parameter estimates and a surrogate measure for the strength of physiological signals (e.g. the BT circadian rhythm). That is, a low uncertainty can be interpreted as corresponding to reliable parameter estimation and a strong physiological signal inherent to the given measurements. In general, the uncertainty estimates were fairly small (table 1; electronic supplementary material, figure S3), suggesting that our method can reliably extract circadian rhythms from wearable data. Importantly, the uncertainty estimate of the phase increased when increasing the number of harmonics (table 1; electronic supplementary material, figure S3). This indicates that the single-harmonic model is more reliable than the multiple-harmonic models in phase extraction. We also found that the population mean of 2.5th percentiles of the amplitude estimates of 24 h rhythm, 12 h rhythm and 8 h rhythm were greater than zero ( p-value < 10 -3 , T-test) (electronic supplementary material, table S4). This indicates that the multiple-harmonic models might be needed to extract all the key physiological information from wearable BT data, matching the previous studies [15,27], while the single harmonic is sometimes enough or more appropriate to identify some parameters, such as the mesor (figure 4f ) and the time of BT minimum (figure 4h and table 1; electronic supplementary material, figure S3).
We next analysed 183 613 days of real-world HR data that were anonymously collected from 1729 subjects by the 'Social Rhythms' Android and iPhone Application (see electronic supplementary material, figure S4) [14,22]. Specifically, we computed the goodness-of-fit statistics and the parameter estimates of the HR data as we did for the BT data (electronic supplementary material, figure S5 and tables S5 and S6). Similar to the result of the BT study (figure 4a-e and table 1), the zero-order harmonic model described the HR data least well. However, there was no significant difference in the goodness-of-fit values between the models (electronic supplementary material, figure S5 and table S5). This modest improvement resulting from consideration of the HR rhythmicity is expected because HR is strongly affected by many non-rhythmic external factors (e.g. stress) [39][40][41] and, compared with them, the intensity of rhythmic signals, for example, from the sinoatrial node of the heart may not be much stronger [14]. We found that the estimates of mesor, phase (defined to be the time of HR minimum), HRpS, and autocorrelation strength did not change significantly when varying the number of harmonics (electronic supplementary material, figure S5 and table S5), and they are similar to the previously reported values [14,42]. The estimates of amplitude increased within the previously reported range (3.96-12.60 bpm) [14,42] when increasing the number of harmonics (electronic supplementary material, figure S5 and table S5), as we observed in the BT data (figure 4 and table 1). Moreover, the uncertainty estimates of the phase increased when increasing the number of harmonics (electronic supplementary material, table S5), suggesting that the single-harmonic model is a suitable choice to reliably extract and analyse the phase.
Using the ALSM with the suitable single-harmonic HR model, we studied the dynamics of the circadian rhythm in HR in free-living conditions. Figure 5 shows the circadian HR phase, defined to be the circadian minimum of HR tracked over a period of greater than 30 days, for four individuals from our Social Rhythms Application dataset. When an individual has consistent daily routines (i.e. regular activity-rest rhythms), the HR phase tracked with sleep (figure 5a). When the activity pattern of individuals suddenly shifts and thus circadian misalignment occurs, the HR phase gradually followed the shift, achieving realignment within 10 days (figure 5b,c). Importantly, we found that there were large interindividual differences in the adjustment to new activity patterns. The HR phase can start to be adjusted to follow the shift in activity pattern without a time delay (figure 5b,c) while it can also remain unchanged (electronic supplementary material, figures S6A and B). Such different patterns can be observed even in an individual (figure 5d). The HR phase remained unchanged at first (figure 5d(i)) and then gradually followed the shift (figure 5d(ii)). These findings matched the previously reported cases obtained by analysing wearable data from ∼ 1000 medical interns [14]. Interestingly, when the activity pattern was advanced right after its delay, increased fluctuations appeared in the HR phase (figure 5d(iii)), showing the circadian disruption effects of chronic external perturbations (e.g. shift work) in realworld settings. Lastly, we sometimes observed a quick shift in the HR phase following a dramatic shift in activity pattern (e.g. around day 50 in figure 5d and electronic supplementary material, figure S7A). These large uncertainties with the quick shifts indicate that accurate parameter estimation becomes challenging when individuals experience continuous external perturbations. This difficulty may arise because evoked components, triggered by external cues, continuously and strongly influence the rhythmic outputs of endogenous oscillators [43], and the influences may not be completely filtered by our method. In this case where oscillatory and evoked components are intricately intertwined, caution should be taken when interpreting the estimation outcomes (electronic supplementary material, figure S7B).
We found that the circadian rhythm in HR typically tracked with sleep (figure 5a-d; electronic supplementary material, figures S6A and B), which cannot be explained by the effects of sleep on HR because all the data collected during sleep were excluded in the analysis as in figure 2a. To further explore the relationship between the circadian rhythm in HR and sleep, we calculated the difference in hours between the estimated HR phase and the sleep midpoint computed using sleep data from wearable devices, denoted by the phase difference. For instance, a phase difference of 2 h means that HR is at its lowest 2 h after the midpoint of sleep. The HR phase closely matched the midpoint of sleep on average (μ = -0.382 h), which was consistent with previous work (μ = -0.45 h) [30] (figure 5e). The range of their differences we found (σ = 2.361 h) was also similar to the previously reported one (σ = 2.25 h) [14]. Moreover, the distribution of average absolute difference that averages only the magnitude of the phase difference (μ = 3.374 h, σ = 1.655 h) closely matched that previously obtained from a different wearable dataset (μ = 3.88 h, σ = 1.56 h) [14] (figure 5f ). We found that the large range of the differences greater than 2 h is due to the inter-and intra-individual difference in the relationship between the HR phase and sleep (figure 5b-d; electronic supplementary material, figures S6C and D). Our findings that are consistent with the previous results [14,22] indicate that circadian assessment can be performed efficiently using our method.
To fully use wearable data when estimating the state space of the molecular clocks in tissues, a data assimilation algorithm based on the KF technique has been recently developed [23]. This integrates the prediction of the clock state simulated by the mathematical model taking wearable activity data as an input with its observation (i.e. the HR phase) (see electronic supplementary material and [23]). It allows us to estimate the evolution of the posterior distributions of the clock phase over days (electronic supplementary material, figure S8). Although this algorithm outperforms the previous one solely based on the phase prediction, its computational cost is high because it extracts the HR phase using the previously developed Bayesian MCMC algorithm (figure 2f ) [14]. We addressed this problem by replacing the Bayesian MCMC method with the computationally efficient ALSM in the KF framework (figure 6a). We tested the updated method on the in silico data, showing that the updated method can successfully estimate the circadian phase (electronic supplementary material, Fig. S8). We next applied it to real-world data collected from a subject showing a large variation in the HR phase (figure 6c) to explore the benefits of the KF approach.
The KF framework with the ALSM allowed us to reduce the uncertainty of the estimate by extracting information from HR data and integrating it with the model prediction (figure 6b-f ). This efficient method provides the possibility of direct analysis of wearable data in consumer-grade wearable devices without resource-intensive steps such as the anonymous transmission of data to secure computer clusters.
Here, we developed a method named ALSM that efficiently extracts physiological information from densely sampled wearable time-series data (figures 1-5). We tested it on over 590 days of BT data and 183 613 days of HR data collected in real-world settings. This allowed systematic comparison of harmonic-regression models, resulting in identifying the suitable models for wearable BT and HR data (figures 3 and 4; electronic supplementary material, figures S3 and S5, and table 1; electronic supplementary material, tables S2-S6). We also showed the usefulness of our method in the assessment of inter-and intra-individual differences in the HR clock (figure 5) and circadian biomarker development (figure 6).
Our work and other ongoing large population field studies such as PiCADo [44], inCASA [45], Intern Health Study [46] and Social Rhythms [14] will rapidly facilitate the validation of wearables.
The extraction of useful interpretable information about physiological processes from wearable data is a longstanding challenge. Approaches for this can be primarily classified into non-parametric and parametric classes. In the non-parametric analysis, wearable data are analysed with non-parametric statistics such as interdaily variability that quantifies the disruption of rhythms [47]. Although new useful nonparametric variables continue to be developed, we are still far from the complete exploitation of wearable data because of challenges in uncertainty quantification of the variables, which is important for clinicians to make a better decision about patient's treatment [48]. To compensate for such limitations, sophisticated parametric methods for wearable physical activity time-series data, such as hidden Markov modelling [48], have recently been proposed. Together with the frameworks for activity data, our Fourier-based method for wearable BT and HR data opens the possibility of increasing the clinical applicability of wearable data. In particular, our method can be used to identify disease-induced abnormal changes in BT and HR by integrating it with machinelearning techniques [30]. The efficiency of our method makes it particularly suited for such integrations.
The HR model (equation (2.5)) directly accounts for the effects of activity on HR, and indirectly accounts for other Here, the analysed real-world wearable-device data were adopted from [23].
Downloaded from https://royalsocietypublishing.org/ on 26 November 2025