You have access

2023Efficient assessment of real-world dynamics of circadian rhythms in heart rate and body temperature from wearable dataJ. R. Soc. Interface.2020230030http://doi.org/10.1098/rsif.2023.0030

Section

Supplemental Material

You have access

Research articles

Efficient assessment of real-world dynamics of circadian rhythms in heart rate and body temperature from wearable data

Dae Wook Kim^†

Dae Wook Kim

http://orcid.org/0000-0003-1060-891X

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

[email protected]

Contribution: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Caleb Mayer^†

Caleb Mayer

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Formal analysis, Funding acquisition, Investigation, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Minki P. Lee

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Investigation, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Sung Won Choi

Division of Pediatric Hematology/Oncology, Department of Pediatrics, University of Michigan, Ann Arbor, MI 48109, USA

Rogel Comprehensive Cancer Center, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Resources, Supervision, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Muneesh Tewari

Rogel Comprehensive Cancer Center, University of Michigan, Ann Arbor, MI 48109, USA

Division of Hematology and Oncology, Department of Internal Medicine, University of Michigan, Ann Arbor, MI 48109, USA

Department of Biomedical Engineering, University of Michigan, Ann Arbor, MI 48109, USA

Center for Computational Medicine and Bioinformatics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Resources, Supervision, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

and

Daniel B. Forger

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

Center for Computational Medicine and Bioinformatics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Dae Wook Kim^†

Dae Wook Kim

http://orcid.org/0000-0003-1060-891X

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

[email protected]

Google Scholar

Find this author on PubMed

Search for more papers by this author

Caleb Mayer^†

Caleb Mayer

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Formal analysis, Funding acquisition, Investigation, Writing – original draft, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Minki P. Lee

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Investigation, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Sung Won Choi

Division of Pediatric Hematology/Oncology, Department of Pediatrics, University of Michigan, Ann Arbor, MI 48109, USA

Rogel Comprehensive Cancer Center, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Resources, Supervision, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Muneesh Tewari

Rogel Comprehensive Cancer Center, University of Michigan, Ann Arbor, MI 48109, USA

Division of Hematology and Oncology, Department of Internal Medicine, University of Michigan, Ann Arbor, MI 48109, USA

Department of Biomedical Engineering, University of Michigan, Ann Arbor, MI 48109, USA

Center for Computational Medicine and Bioinformatics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Resources, Supervision, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

and

Daniel B. Forger

Department of Mathematics, University of Michigan, Ann Arbor, MI 48109, USA

Center for Computational Medicine and Bioinformatics, University of Michigan, Ann Arbor, MI 48109, USA

Contribution: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review & editing

Google Scholar

Find this author on PubMed

Search for more papers by this author

Published:23 August 2023https://doi.org/10.1098/rsif.2023.0030

Abstract

Laboratory studies have made unprecedented progress in understanding circadian physiology. Quantifying circadian rhythms outside of laboratory settings is necessary to translate these findings into real-world clinical practice. Wearables have been considered promising way to measure these rhythms. However, their limited validation remains an open problem. One major barrier to implementing large-scale validation studies is the lack of reliable and efficient methods for circadian assessment from wearable data. Here, we propose an approximation-based least-squares method to extract underlying circadian rhythms from wearable measurements. Its computational cost is ∼ 300-fold lower than that of previous work, enabling its implementation in smartphones with low computing power. We test it on two large-scale real-world wearable datasets: $\sim 600 days$ of body temperature data from cancer patients and ∼ 184 000 days of heart rate and activity data collected from the ‘Social Rhythms’ mobile application. This shows successful extraction of real-world dynamics of circadian rhythms. We also identify a reasonable harmonic model to analyse wearable data. Lastly, we show our method has broad applicability in circadian studies by embedding it into a Kalman filter that infers the state space of the molecular clocks in tissues. Our approach facilitates the translation of scientific advances in circadian fields into actual improvements in health.

1. Introduction

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⁻¹) 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.

Figure 1. The approximation-based least-squares method (ALSM) can estimate parameters that account for the influence of circadian and state-evoked external components on physiological outputs from wearable data. (a) Wearables allow continuous measurement of physiological signals such as body temperature (BT), heart rate (HR) and activity with high-frequency resolution. (b) The ALSM can extract the circadian rhythms in BT and HR and the effect of activity on HR from *in silico* wearable data of BT (upper panel), HR and activity (lower panel) measured at 1 min intervals. The black solid line and dashed lines in (b) represent the mean, 2.5% and 97.5% quantiles of the estimated time traces. (c,d) The samples of parameter estimates obtained with the ALSM using the sample time trace of BT (c) and HR (d) shown in (b). 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² 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].

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. 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 activity-derived 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 $\sim 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–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-plus-first-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).

y_{t} = s_{t} + v_{t},

2.1

where

\begin{aligned} s_{t} & = μ + \sum_{r = 1}^{n} [a_{r} \cos (\frac{2 π r}{τ} t) + b_{r} \sin (\frac{2 π r}{τ} t)], \\ v_{t} & = α \cdot v_{t - Δ t} + ϵ_{t}, \end{aligned}

and it is assumed that the circadian period τ = 24 h, |α| < 1, and the

ϵ_{t}

values are distributed as Gaussian random variables with mean zero and variance

σ_{ϵ}^{2}

. The autoregressive noise process describes the ongoing effects of external factors on the BT. The Gaussian noise

ϵ_{t}

represents new external influences, for example, from stress, hormones and measurement error. The autocorrelation factor α represents the ongoing contribution of external factors.

2.1.2. Approximation-based linear least-squares method to estimate the body temperature model parameters

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:

\begin{aligned} y_{t} & = s_{t} + v_{t} \\ = s_{t} + α \cdot v_{t - Δ t} + ϵ_{t} \\ = (1 - α) \cdot s_{t} + α \cdot (s_{t} + v_{t - Δ t}) + ϵ_{t} \\ \approx (1 - α) \cdot s_{t} + α \cdot (s_{t - Δ t} + v_{t - Δ t}) + ϵ_{t} \\ = (1 - α) \cdot s_{t} + α \cdot y_{t - Δ t} + ϵ_{t} . \end{aligned}

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)^′, $θ = {(\tilde{μ}, {\tilde{a}}_{1}, {\tilde{b}}_{1}, \dots, {\tilde{a}}_{n}, {\tilde{b}}_{n}, α)}^{'}$ and $ϵ = {(ϵ_{Δ t}, ϵ_{2 Δ t}, \dots, ϵ_{t - Δ t}, ϵ_{t})}^{'}$ where $\tilde{μ} = (1 - α) \cdot μ$ , ${\tilde{a}}_{r} = (1 - α) \cdot a_{r}$ , and ${\tilde{b}}_{r} = (1 - α) \cdot b_{r}$ . Then, equation (2.2) can be rewritten in matrix notation as

y = X θ + ϵ,

2.3

where

X = [\begin{matrix} 1 & \cos (\frac{2 π Δ t}{τ}) & \sin (\frac{2 π Δ t}{τ}) & \dots & \cos (\frac{2 π n Δ t}{τ}) & \sin (\frac{2 π n Δ t}{τ}) & y_{0} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋮ \\ 1 & \cos (\frac{2 π t}{τ}) & \sin (\frac{2 π t}{τ}) & \dots & \cos (\frac{2 π n t}{τ}) & \sin (\frac{2 π n t}{τ}) & y_{t - Δ t} \end{matrix}] .

By applying the linear LSM to equation (2.3), we can derive the probability density function for

\hat{θ}

that minimizes ||y − Xθ|| given y and X, where || · || denotes the 2-norm (equation (2.4)).

\hat{θ} \sim N ({(X^{'} \cdot X)}^{- 1} \cdot X^{'} \cdot y, {(X^{'} \cdot X)}^{- 1} \cdot η^{2}),

2.4

where

η^{2} = \frac{1}{dim (y) - dim (θ) - 1} \cdot (Γ^{'} \cdot Γ)

and

Γ = y - X \cdot {(X^{'} \cdot X)}^{- 1} \cdot X^{'} \cdot y,

and

dim (\cdot)

represents the dimension of a vector. Then, we can compute the probability density of parameters μ, a_r , b_r and α from the distribution of

\hat{θ}

(equation (2.4)) using the equations

μ_{r} = {\tilde{μ}}_{r} / (1 - α)

a_{r} = {\tilde{a}}_{r} / (1 - α)

and

b_{r} = {\tilde{b}}_{r} / (1 - α)

. From these, we can compute the mean and uncertainty estimates of the phase and amplitude of the signal s_t (electronic supplementary material).

2.2. Heart rate parameter estimation from wearable data

2.2.1. A harmonic-regression model of the human heart rate rhythm

The human HR shows a circadian variation [17]. This HR circadian signal has also been successfully analysed with a harmonic-regression-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:

y_{t} = s_{t} + d \cdot a_{t} + v_{t},

2.5

where a_t is the activity level (i.e. step count) at time t and d is the increase in HR per step (HRpS).

2.2.2. Approximation-based nonlinear least-squares method to estimate the heart rate model parameters

We next describe the nonlinear version of the approximation-based LSM (ALSM) to efficiently estimate the HR model parameters. Like the linear version, we first reformulate and approximate equation (2.5) as follows:

\begin{aligned} y_{t} & = s_{t} + d \cdot a_{t} + v_{t} \\ = s_{t} + d \cdot a_{t} + α \cdot v_{t - Δ t} + ϵ_{t} \\ = (1 - α) \cdot s_{t} + α \cdot (s_{t} + d \cdot a_{t - Δ t} + v_{t - Δ t}) + d \cdot a_{t} \\ - α d \cdot a_{t - Δ t} + ϵ_{t} \\ \approx (1 - α) \cdot s_{t} + α \cdot (s_{t - Δ t} + d \cdot a_{t - Δ t} + v_{t - Δ t}) + d \cdot a_{t} \\ - α d \cdot a_{t - Δ t} + ϵ_{t} \\ = (1 - α) \cdot s_{t} + α \cdot y_{t - Δ t} + d \cdot a_{t} - α d \cdot a_{t - Δ t} + ϵ_{t} . \end{aligned}

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 $θ = {(\tilde{μ}, {\tilde{a}}_{1}, \dots, {\tilde{a}}_{n}, {\tilde{b}}_{1}, \dots, {\tilde{b}}_{n}, \tilde{d}, α)}^{'}$ , denoted by $\bar{θ}$ , and its uncertainty (i.e. the covariance matrix $Σ_{θ}$ ) by using any LSMs for nonlinear squares fitting problems. In this study, we used the Levenberg–Marquardt algorithm (see electronic supplementary material and [31]).

3. Results

3.1. Extraction of physiological parameters from wearable data

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 $\sqrt{a_{r}^{2} + b_{r}^{2}} \sin (π r / 24) Δ t$ 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].

Figure 2. The ALSM can track the circadian HR phase much more efficiently than the previous method although there is a misalignment between the HR and activity patterns. (a,b) Although the HR data collected during the wake period are only used for estimation (a), the ALSM can successfully infer the physiological parameters (b). The estimation was done as in figure 1b,d. (c) Although the measurement interval is increased, the estimation is still reasonably successful. We fitted the HR model to 2-day intervals cent red at the sleep time in between, as in (a). 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.

3.2. Efficient tracking of the circadian heart rate phase from wearable data

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_{i}^{off}$ and onset time $t_{i}^{on}$ of the day i were defined to be $t_{i}^{off} \sim 7.00 + N (0, σ_{t}^{2})$ and $t_{i}^{on} \sim 23.00 + N (0, σ_{t}^{2})$ 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 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.

3.3. The circadian phase can be extracted from 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³ 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 single-harmonic 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.

3.4. Extraction of circadian physiological parameters from body temperature data

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-mean-square error (RMSE), the coefficient of determination denoted by R², 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² 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.

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.
Collapse
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²	−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 ( $^{\circ} C$ )	35.983 ± 0.639	35.985 ± 0.636	35.985 ± 0.636	35.985 ± 0.636
amp. ( $^{\circ} 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 ( $^{\circ} C$ )	0.240 ± 0.142	0.168 ± 0.092	0.142 ± 0.078	0.124 ± 0.067
amp. ( $^{\circ} 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

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 two-harmonic 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 two-harmonic 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^{\circ}$ C and $37^{\circ} C$ , which closely matched the normal range of the temperature of the skin surface of the trunk ( ${33.5}^{\circ} C - {36.9}^{\circ} 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}^{\circ} C$ ) was consistent with the previously reported one ( ${0.88}^{\circ} 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 $\sim {0.6}^{\circ} C$ regardless of the number of harmonics used for estimation (electronic supplementary material, table S4), which closely matched the previously reported value (approx. ${0.5}^{\circ} 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⁻³, 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).

3.5. Estimation of real-world dynamics of the circadian rhythm in heart rate

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–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 real-world 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.

3.6. Efficient estimation of the clock state

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.

Figure 6. Estimation of the circadian phase of the molecular clocks in free-living conditions. (a) Diagram illustrating the KF method for the clock phase estimation. It integrates the model prediction from wearable activity data (i) with the peripheral HR phase estimated using our ALSM (ii). As a result, the clock phase predicted by the limit-cycle model [21] is adjusted, which typically improves the estimation accuracy (electronic supplementary material, figure S8). (b–d) The clock phase estimates. The clock phase can be predicted by the limit-cycle model taking activity as an input (b). This prediction can be updated with the HR phase estimated by the ALSM (c) using the KF algorithm (a), resulting in the corrected clock estimate (d). The graphs were plotted as in figure 5a. (e) The probability density distributions of the clock phase estimated by the three methods on the last day. (f) The phase uncertainty defined by the standard deviation of the phase estimate over days. Here, the analysed real-world wearable-device data were adopted from [23].

4. Conclusion

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 non-parametric 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 machine-learning 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 changes to HR on shorter timescales through the autocorrelation. Previous work has shown the success of this model in distinguishing HR phase estimates from sleep- or activity-derived timing metrics, and generating parameter estimates that match those from a constant routine [14,22]. The parameters estimated from our new algorithm are also comparable to those from this previous work. Furthermore, the impact of posture on HR has been tested using a similar modelling framework [30], and the derived signal has been shown to act similarly to an internal clock [23]. These findings suggest that the phase and amplitude parameters extracted from the model correspond to real aspects of circadian timekeeping. However, the possibility remains for real-world confounders, particularly as this method becomes applied to a wider range of datasets and populations. Therefore, large-scale validation of this algorithm against a gold-standard circadian marker or other meaningful external data source could be valuable to confirm these results in the future. We showed that our method produced comparable results to the previous method (figure 5), which has shown consistent performance across various devices such as the Apple Watch, Fitbit and Mi Band [14]. This provides evidence supporting the consistency of our method across devices. It would be interesting to explore the performance of our method across different wearables by comparing estimates obtained from individuals wearing two devices simultaneously, as done in previous work [14].

When extracting the circadian phase from noisy HR data with gaps, the single harmonic model generally outperforms the multiple harmonic models (figure 3). Counterintuitively, this remains the case even when the underlying data is generated using multiple harmonics. We note that while including additional harmonic terms in the model will reduce the mean absolute error between the fit and the data points, it does not guarantee an improvement in the phase extraction, as the phase error itself is not minimized. Indeed, the true phase in wearable HR data is often obscured by noise and gaps in the measurements. Hence including higher-order harmonic terms in the model may allow overfitting of the noise, particularly under the assumption that the higher-order harmonic rhythms have a relatively small amplitude compared with the noise and 24 h hour rhythm. Indeed, the single harmonic model outperforms the higher-order models most clearly under challenging data conditions (i.e. large gap lengths and measurement intervals; electronic supplementary material, table S2). Thus, the phase prediction may suffer when the model cannot distinguish between the contributions of the higher-order harmonics and noise components in the presence of gaps. Therefore, the simpler model that considers the dominant 24 h rhythm yields sufficient or even optimal performance in many cases.

In this study, we analysed axillary temperature from subjects who were hospitalized, providing a partially controlled environment (e.g. physical activities were limited), using a wearable thermometer, TempTraq. Recent work showed that TempTraq axillary temperature measurements were highly correlated to measurements of tympanic temperature [49], which is considered a core temperature [50]. Moreover, both classic and recent studies support the accuracy of axillary temperature measurements as reflections of core body temperature measurements [51,52]. This indicates that the non-invasive collection of high-frequency axillary temperature measurements using TempTraq from hospitalized subjects offers a promising approach for studying the body temperature rhythm.

Currently, wearable measurements are transmitted to external computer clusters for circadian assessment with the Bayesian MCMC framework, which is computationally demanding [14]. Then, the circadian parameters extracted from the data are sent back to wearable devices. These intermediate steps are not only labour and resource intensive but also could lead to data breach issues [53]. One way to address this is to directly perform the assessment only using the built-in computing power of wearables and smartphones. For this, algorithms that are significantly more computationally efficient than the previous ones need to be developed. In this circumstance, our computationally fast approach (figure 2f) provides a promising starting point.

Here, we proposed a framework allowing efficient estimation of parameters of any regression models with correlated noise. We used this approach to analyse wearable data by adopting the harmonic-regression-plus-first-order-autoregressive models as in the previous studies [14,15,22,30]. Brown and Czeisler showed that the first-order autoregressive noise process is enough to analyse core BT data collected in laboratory settings [15], which matches our findings (figures 3 and 4; electronic supplementary material, figures S3 and S5, and table 1; electronic supplementary material, table S5). It would be interesting to investigate whether this conclusion is preserved across different wearable datasets using other methods [14,15] as well as the ALSM that can be extended to harmonic models with high-order correlated noise (electronic supplementary material).

Data accessibility

The computer codes used in this study are accessible from the GitHub repository: https://github.com/daewookk/ALSM. We adopted the temperature wearable data from [16] by obtaining permission from the authors of the previous study. The raw heart rate and activity; wearable data were collected under the terms of Social Rhythms Privacy Policy, as done in the previous work [22], which allows the sharing of summary data. The summary data of individual subjects presented in the study (electronic supplementary material, figure S4) are deposited from the GitHub repository: https://github.com/daewookk/Socialrhythms_J_R_SOC_Inter.

The data are provided in electronic supplementary material [54].

Declaration of AI use

We have not used AI-assisted technologies in creating this article.

Authors' contributions

D.W.K.: conceptualization, data curation, formal analysis, investigation, methodology, project administration, resources, software, supervision, validation, visualization, writing—original draft, writing—review and editing; C.M.: formal analysis, funding acquisition, investigation, writing—original draft, writing—review and editing; M.P.L.: investigation, writing—review and editing; S.W.C.: resources, supervision, writing—review and editing; M.T.: resources, supervision, writing—review and editing; D.B.F.: conceptualization, funding acquisition, project administration, supervision, writing—review and editing.

All authors gave final approval for publication and agreed to be held accountable for the work performed therein.

Conflict of interest declaration

D.B.F. is the CSO of Arcascope, a company that makes circadian rhythms software. Both he and the University of Michigan own equality in Arcascope.

Funding

This work was funded by Human Frontiers Science Program Organization grant no. RGP0019/2018, NSF DMS grant no. 2052499 and ARO MURI grant no. W911NF-22-1-0223. C.M. is grateful for support through the SBIR grant no. R43CA236557. S.W.C. and M.T. acknowledge support from A. Alfred Taubman Medical Research Institute Grand Challenge and Innovation Project grants. SWC also acknowledge support from NIH/NHLBI grants R01HL146354 and K24HL156896 and NCI grant R01CA249211.

Acknowledgements

We thank the LSA Technology Services team of the University of Michigan for their technical support.

Footnotes

^†These authors contributed equally to the study.

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.6778066.

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Journal of The Royal Society Interface cover image

August 2023

Volume 20Issue 205

Article Information

DOI:https://doi.org/10.1098/rsif.2023.0030
PubMed:37608712
Published by:Royal Society
Online ISSN:1742-5662

History:

Manuscript received23/01/2023
Manuscript accepted31/07/2023
Published online23/08/2023

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

Efficient assessment of real-world dynamics of circadian rhythms in heart rate and body temperature from wearable data

Abstract

1. Introduction

2. Methods

2.1. Body temperature parameter estimation from wearable data

2.1.1. A harmonic-regression model of the human body temperature rhythm

2.1.2. Approximation-based linear least-squares method to estimate the body temperature model parameters

2.2. Heart rate parameter estimation from wearable data

2.2.1. A harmonic-regression model of the human heart rate rhythm

2.2.2. Approximation-based nonlinear least-squares method to estimate the heart rate model parameters

3. Results

3.1. Extraction of physiological parameters from wearable data

3.2. Efficient tracking of the circadian heart rate phase from wearable data

3.3. The circadian phase can be extracted from wearable data more accurately with a single-harmonic model than with a multiple-harmonic model

3.4. Extraction of circadian physiological parameters from body temperature data

3.5. Estimation of real-world dynamics of the circadian rhythm in heart rate

3.6. Efficient estimation of the clock state

4. Conclusion

Data accessibility

Declaration of AI use

Authors' contributions

Conflict of interest declaration

Funding

Acknowledgements

Footnotes

This Issue