WO2020074881A1 - Resting heart rate estimation - Google Patents

Abstract

A method, computer program and device for determining a user's resting heart rate are disclosed. The method comprises: receiving data indicative of a user's heart rate sensed over a period of time; fitting at least some of the sensed data to a base pattern indicative of estimated variations in an average person's resting heart rate over a predetermined time period, the base pattern comprising a cosinor comprising a plurality of harmonics; and determining the resting heart rate from the sensed data fitted to the base pattern.

Description

RESTING HEART RATE ESTIMATION

FIELD OF THE INVENTION

The field of the invention relates to the estimation of a resting heart rate. BACKGROUND

There are an increasing number of devices being used to monitor things such as an individual’s heartbeat. Various information can be derived from these measurements that may be of interest to a user. For example, changes in heart rate can be used to determine various things such as how hard a user is working when exercising and the number of calories burnt. However, different users have different physiology and thus, such calculations should be tailored to the user. One factor that is important in these estimations is the user’s usual resting heart beat, however, it is not straightforward to accurately measure this. A user’s resting heart rate will vary depending on the time of day and the activity levels of the user and thus, determining an average resting heart rate that accurately characterises a user can be difficult.

It would be desirable to be able to determine a user’s resting heart rate in a

computationally efficient and accurate manner.

SUMMARY

A first aspect provides a method of determining a user’s resting heart rate, comprising: receiving data indicative of a user’s heart rate sensed over a period of time; fitting at least some of said sensed data to a base pattern indicative of estimated variations in an average person’s resting heart rate over a predetermined time period, said base pattern comprising a cosinor comprising a plurality of harmonics; and determining said resting heart rate from said sensed data fitted to said base pattern.

There are an increasing number of devices available to continuously or quasi continuously measure a user’s heart rate and this information provides insights not only into how strenuous an activity may be for a user but also on their current wellbeing and their physical fitness. A user’s heart rate is affected by many things such as activity level, mood, stress, eating and sleep. A user’s resting heart rate provides information regarding their fitness and also serves as a base line against which any increase in heart rate due to other factors can be measured. However, as a user’s heart rate varies throughout the day, determining a resting heart rate is not straightforward and may be mathematically complex and computationally expensive. The first aspect of the present invention addresses these issues by fitting data indicative of a user’s heart rate to a base pattern that provides a pattern indicative of estimated variations in an average person’s resting heart rate. This base pattern is formed from a cosinor with a plurality of harmonics. As noted previously a person’s resting heart rate will vary over a period of time and this variation may be periodic to some degree in its nature, as users sleep, work and eat at similarly spaced points in a day. Thus, modelling the variations in the resting heart rate using a cosinor with multiple harmonics and fitting measured heart rate to a base pattern modelled in this way, is mathematically efficient while providing a useful representation of how heart rate generally varies. Thus, a method of determining a resting heart rate that is efficient in processing power and provides an effective output is provided.

In this regard the resting heart rate is the heart rate where data while active is filtered out. The actual heart rate is the variable that is measured and the resting heart rate is estimated from this.

In some embodiments, the method further comprises sensing the heart beat and receiving the heart beat data from the sensor, while in other embodiments the step of receiving comprises receiving the data from a remote sensor.

Although the predetermined time period over which the measurements are analysed may be a number of things, in some embodiments said predetermined time period comprises 24 hours. People have routines that vary over a 24 hour period and thus, this is generally the most effective time period to analyse a user’s heart beat variations over.

Although the plurality of harmonics may comprise two or three harmonics in some embodiments said plurality of harmonics comprise four harmonics, said four harmonics comprising one day, a half day, a quarter of a day and an eighth of a day.

It has been found that four harmonics over a 24 hour period provides a cosinor that is particularly effective for mapping the resting heart rate of a person over this period.

In some embodiments, said sampling is performed over several days and said method further comprises the step of combining heart rate data for a same period of time from different days in said sampling period. Although, the data may be sensed over a single predetermined time period such as 24 hours, in some embodiments it is sensed over several days and the method comprises the step of combining heart rate data for a same period of time from different days in said sampling period. As the user’s heart rate is often continuously or quasi- continuously sampled, the data from several time periods may be combined to improve the accuracy of the estimation.

In some embodiments, said step of combining comprises weighting said data prior to said combining, higher heart rate data being given a lower weighting.

As it is the resting heart rate that is of particular interest then it may be advantageous if lower heart rate data is given a higher weighting than higher heart rate data. Higher heart rate data may be indicative that the user is not currently resting but is involved in some activity and thus, the lower heart rate data is probably more typical of a resting heart rate.

In some embodiments, the method further comprises determining at least one of sex and age of said user and adjusting said base pattern in dependence upon said at least one of said sex and age of said user prior to performing said fitting step.

There are some factors which generally affect the resting heart rate of a person and where these are known they can be used to modify the base pattern to form a base pattern that is more likely to be characteristic of a user of that age or sex. In this way the base pattern is improved for that user and the fitting step becomes easier and more accurate.

In some embodiments, the method further comprises determining said user’s activity levels and selecting heart rate data during periods of time where activity level is low and discarding heart rate data during periods of time where activity level is deemed to be high.

As it is the resting heart rate that is being determined discarding heart rate data when a user’s activity levels are deemed to be increased renders the sampled data more characteristic of the heart rate being determined and improves the accuracy of the estimation. Users’ activity levels may be determined in a number of ways, in some cases there may be a sensor for sensing a user’s activity levels such as a step detector whose input is received alongside the heart rate data and can be used to discard data where the activity levels are high. In other cases, changes in heart rate maybe monitored and maybe used to determine where they rise that an activity has probably commenced and these higher levels can be discarded.

In some embodiments, said method comprises applying a weighting to at least some of said user data prior to performing said combining step.

Weighting can be applied to the user data to assign more importance to some of said data than to others.

The model used to estimate resting heart rate is such that each data point used in the fitting step may be assigned a weight prior to the fitting step. This allows a more advanced analysis to be performed where the importance of certain data points can be uplifted or downgraded depending on circumstances. For example, a higher heart rate is indicative of a person moving and as such is less characteristic of the resting heart rate than other data. It may therefore be advantageous to apply a higher weight to lower heart rates than to higher heart rates. Furthermore, where for example, a real time streaming analysis is being performed then this may be achieved by appropriate weighting of data so that the method may apply a higher weight to newer more recently measured data than to older data. This allows the model to continually update the determination, and use of selected weighting values can be used to control this process. In some embodiments, said step of fitting comprises forcing said user’s data to have a same basic shape as said base pattern, said same basic shape having a freedom of a vertical translation, a scale factor and a horizontal translation with respect to said base pattern.

A computationally efficient way of fitting the data is by forcing the user’s data to have the same basic shape as the base pattern but providing it with the freedom of vertical translation, a scale factor and a horizontal translation with respect to the base pattern. Not only is this computationally efficient but it also provides an accurate indication of the user’s rest heart rate based on the assumption that users have similar periodic variations in their heart rate but that the resting heart rate may be a factor higher for some people, they may have larger differences between day and night resting rates and their sleep patterns may occur at different times of the day. These can be accounted for by the freedoms provided in fitting. In effect the refined user model is forced, after fitting the 3 parameters that are known to change, to have the same "shape" as the expected shape of circadian effect on heart rate. Thus, when fitting the model the only parameters that are allowed to change are those that we know from medical literature we can change. Average HR is linked to the general fitness of the user, amplitude is again linked to the general fitness of the user, phase is linked to when you go to sleep and wake up. We don't allow the model to change "shape" because it's known that circadian effect on HR has a certain shape, so that shape is forced. This renders the calculation to determine the resting heart rate both simpler and more accurate.

In some embodiments, said step of fitting comprises forcing said user’s data to have a same basic shape as said base pattern by performing harmonic regression on said data to force the shape of the variation in heart rate over time to be substantially the same as the base pattern, any variations between the base pattern and user’s pattern being in the scale factor of the shape the vertical translation of the shape , and in the horizontal translation or phase offset of the shape with respect to the base pattern.

Having determined that the shape should be fitted in this manner with these freedoms, harmonic regression can be performed on the data which will force the shape to have substantially the same shape as the base pattern.

In some embodiments, said step of fitting comprises applying user data to the cosinor of the base pattern and determining a cost function by squaring the difference between measured heart rate and modelled heat rate weighted to favour lower heart rates, and determining an analytic solution that minimises the cost function.

When performing the fitting, one way of doing so that is both accurate and efficient is to determine the modelled heart rate for that user and compare it with the measured heart rate and square the difference to find the lowest difference which would be the most accurately modelled solution. In addition this difference can be weighted to favour lower heart rates an analytic solution then being found that minimises the cost function.

In some embodiments, said cost function to be minimised comprises a further regularisation term that adds a cost where the cosinor with the user’s data is far from the base pattern. Additionally further regularisation can be applied to the cost function this

regularisation adding a cost where the cosinor formed with the user data is far from the base pattern. This favourises data that is close to that of the model.

In some embodiments, the vertical translation of the basic shape to the base pattern provides an indication of amplitude difference in heart rate between night and day, the scale factor provides an indication of average resting heart rate and the horizontal translation provides an indication of sleep times of the user.

As noted previously there are various freedoms provided when fitting the shape to the base pattern and these correspond to various characteristics of the user’s heart rate. Thus, the amplitude difference is related to the difference in heart rate between night and day of a user compared to an average user, while the scale factor provides an indication of the average resting heart rate and the horizontal translation provides an indication of sleep times of the user and how they differ from the model which corresponds to the average user.

In some embodiments, said step of determining said resting heart comprises determining at least one of an average resting heart rate and a resting heart rate as a function of time of day.

Although, the resting heart rate can be determined as an average resting heart rate from the scale factor of the model for example, as explained above, the resting heart rate can also be determined as a function of time of day from the model itself. Both of these can be of interest to the user, the average resting heart rate providing an indication of fitness of a user for example, while the resting heart rate for the time of day being useful when trying to determine the activity of a user for example and how much their heart rate might have been increased by that activity.

In some embodiments, the method comprises a further step of outputting an indication of said resting heart rate, wherein said indication comprises at least one of an average value of a determined resting heart rate, a current resting heart rate and an indication of fitness.

In some embodiments, said at least some of said data comprises two surrogates derived from a measure indicative of a current resting heart rate and a time of the measure. The method provides a computation that does not require a full data set but can perform the resting heart rate determination based on only two surrogates, these two surrogates being values derived from a measure of the resting heart rate and the time of the measure. This means the method does not require access to the full dataset, but only a limited set of values. All these values are moreover an additive function of data sets, which implies that they can be computed separately for each subset of samples, and then added together to obtain the total value. This allows the data to be accumulated directly in the user’s device, with big advantages in computational complexity, data transfer and privacy. The model is computationally cheap, depending linearly on the size of the data.

Furthermore, the method is such that the fitting of the data goes beyond the fitting of circadian activity to a resting heart rate and can be used not just to fit real, floating point data to a pattern, but to fit vectorial data such as complex data to different periodic patterns. The vectorial data may include data from an accelerometer or gyroscope in addition to the heart rate data and can be processed by this method to fit data to non-periodic patterns using a linear combination of an arbitrary set of functions, for example, allowing a determination to be performed in multiple dimensions for multiple features.

In some embodiments, said two surrogates comprise a shape coefficient S and Fourier convolution coefficient F.

The computation of the model does not need the full dataset, but in some

embodiments, only two surrogates (F and S), derived where the data is accumulated. This implies that the model can be implemented in a streaming approach, where the data is accumulated on the device, with important consequences for security and privacy of the data that never leaves the user devices.

In some embodiments, the Fourier convolution coefficient F_d = åjWjyje^i2Itdx-i and the shape coefficient S = åjWje^i2,tdx-i. Where y, is the j^th measure of RHR of the user and Xj is the time of the measure and W_j is the weight applied to the j^th measure.

In some embodiments, said at least some of said data is stored on said device.

As the determination can be performed on a subset of the full data set, the data accumulated on the device may be a corresponding subset allowing the other data to be discarded, this allows the model to be implemented in a streaming approach where the reduced data is accumulated on the device. This has important consequences for security and privacy of the data as it need never leave the device.

A second aspect provides a method of comparing a user’s heart rate to a predicted value, said method comprising: sensing a user’s activity level and determining an expected heart rate from an estimated increase in heart rate due to said sensed activity level and a determined resting heart rate for that time of day, said resting heart rate for said time of day being determined by a method according to a first aspect; outputting an indication of a difference in measured heart rate and predicted heart rate.

As noted above determining the resting heart rate for a period of time allows a user to determine from a current heart rate how much it has increased, the increase can then be used to determine how hard the user is working and perhaps as an indication of calories burnt.

In some embodiments, the method comprises determining from said difference in measured heart rate and predicted heart rate a likelihood of said measured heart rate occurring.

In addition to determining the difference in the heart rate from the predicted heart rate, the probability or likelihood of that heart rate for a user can be determined, the closer the value being to the predicted or expected value, the more likely it is to occur,

In some embodiments, the method comprises determining a plurality of said likelihoods of a plurality of measured heart rates occurring, said plurality of measured heart rates being measured over a predetermined time period, and determining an average of said plurality of determined likelihoods during said time period.

In some embodiments, where said averaged difference is greater than a predetermined value, said indication comprises an alert.

Determining an average of likelihoods of heart rate occurring over a time window, of anything between 30 minutes and 6 hours, has been found to be a particularly accurate way of detecting when a heart rate is“very strange”. This can be used to trigger an alert and can indicate interesting things such as stressful days, stressful meetings, illness, sleep disruptions etc. Although the average can calculated in a number of ways in some embodiments, the average comprises a geometric mean of said likelihood values for said time period.

In some embodiments, said indication comprises a wellbeing indication.

The current resting heart rate may also be an indication of wellbeing an increased resting heart rate perhaps indicating either anxiety or an underlying illness.

A third aspect of the invention provides a method of determining a base pattern indicative of resting heart rate variations in an average person over a period of time, comprising: receiving data from a plurality of users indicating heart rate over said period of time; selecting heart rate data during periods of time where activity level is deemed to be low and discarding heart rate data during periods of time where activity level is deemed to be high; combining said data from multiple users and multiple days; generating a base pattern indicative of estimated changes in resting heart rate over time for an average person by fitting a periodic model to said combined data, said periodic model comprising a cosinor with a plurality of harmonics.

In order to accurately determine a user’s resting heart rate a base pattern is used against which the user’s heart variations over a time period are compared. In order for the model to work well, the base pattern should have shape that is indicative of a user’s heart rate variations over time. One way of providing an effective and accurate base pattern is to generate it from data received from a plurality of users that is filtered by selecting a heart rate during periods of time when activity level is deemed to be low and combining this data for the multiple users and the multiple days and from this data generating the base pattern. The base pattern is formed from a periodic model which comprises cosinor with a plurality of harmonics. In this regard the inventors recognise that there is a periodic nature to heart rate variations and as such modelling it with a cosinor with a plurality of harmonics allows it be represented in a relatively accurate and yet computationally efficient manner.

In some embodiments said step of combining comprises an initial step of averaging heart rate data from each user for a same period of time in a day from a sample period of several days and coalescing said data for a plurality of users. Where the heart rate data being collected over several days then for each user a heart rate pattern for that user is determined by averaging the data for those days and then the data for the plurality of users is coalesced.

In some embodiments, said step of coalescing comprises combining said data, said data being weighted prior to said combination, higher heart rate data having a lower weighting.

Activities and external factors can increase the heart rate of a user and thus, in order to improve accuracy and make the data more representative of the resting heart rate it may be advantageous to weight higher heart rate data so it has a decreased effect on the base pattern. In this regard, where the data from one user over several days is combined over this time period such that an average heart rate of that user for one day is provided then the weighting may be done at this point so that for the days where the heart rate data is higher at a particular period this heart rate data is weighted to provide less of a contribution to the final result than any lower heart rate data sampled for the same period. Additionally and/or alternatively the heart rate data with higher values may be given the lower weighting at the coalescing step where different users’ data is combined. In this regard, heart rate measures are influenced by factors such as mental stress, flu etc. and by providing this weighting more importance is given to points with a low heart rate, which are most likely to be unaffected by these external factors.

In some embodiments, said step of fitting comprises using harmonic regression to said sampled data to fit said sampled data to said shape defined by said cosinor.

In some embodiments, said plurality of harmonics comprise four harmonics, said four harmonics comprising one day, a half day, a quarter of a day and an eighth of a day.

A fourth aspect provides a method of determining a user’s resting heart rate according to a first aspect wherein said base pattern is determined according to a method according to a third aspect.

A fifth aspect provides a computer program comprising computer executable instructions which when executed by a processor are operable to control said processor to perform a method according to a first aspect, a second aspect or a third aspect.

A sixth aspect provides a device for determining a user’s resting heart rate, comprising: a sensor configured to continually sense said user’s heart rate over a period of time; processing circuitry configured to fit at least some of said sensed data to a base pattern indicative of estimated variations in an average person’s resting heart rate over a predetermined time period, said base pattern being modelled by a cosinor comprising a plurality of harmonics and to determine said resting heart rate from said sensed data fitted to said base pattern.

In some embodiments said processing circuitry is configured to estimate values for at least some of said intermittently received data during periods of time that said data is not received and to analyse and combine both said received data and estimated values.

Some of the data received from the sensors may be received continuously or almost continuously while other data is received intermittently. The mathematical tools we developed adapt the received heart rate data of a user and in particular, that deemed to be resting heart rate data to the expected shape taking account of the circadian rhythm effect on resting heart rate.

For example in some cases the biometric data collected may be heart rate data and this may be fitted to a certain base pattern. Usually the data needed for this kind of analysis is sampled at a constant rate, i.e. evenly timed samples, e.g. every 5 minutes, or every 10 minutes. However, our data contains missing datapoints and is unevenly sampled, so we developed the mathematical tools outlined in the description below to enable us to fit data that is not evenly sampled.

Further particular and preferred aspects are set out in the accompanying independent and dependent claims. Features of the dependent claims may be combined with features of the independent claims as appropriate, and in combinations other than those explicitly set out in the claims.

Where an apparatus feature is described as being operable to provide a function, it will be appreciated that this includes an apparatus feature which provides that function or which is adapted or configured to provide that function.

BRIEF DESCRIPTION OF THE DRAWINGS

Embodiments of the present invention will now be described further, with reference to the accompanying drawings, in which: Figure l shows a flow diagram schematically illustrating a method of generating the base pattern;

Figure 2 shows a flow diagram schematically showing a method of determining a user’s resting heart rate;

Figures 3a to 3b shows the distribution of number of HR samples per user and of users per age;

Figure 4 shows HR used for training the prior and the resulting prior for the circadian rhythm;

Figure 5: HR data of user 04E93AE4B9 with prior and posterior;

Figure 6: Four users’ HR data with prior and posterior;

Figure 7: Distribution of circadian parameters;

Figure 8: Number of datapoint (logarithmic scale) by combination of minutes active in the last 5 minutes and in the last 20 minutes;

Figure 9: Residual mean by minutes active in the last 5 minutes and in the last 20 minutes;

Figure 10: Residual standard deviation by minutes active in the last 5 minutes and in the last 20 minutes;

Figure 11: Residual means estimated with a linear function;

Figure 12: Distribution of activity parameters;

Figure 13: Distribution HR prediction relative error;

Figure 14: HR of user 11 with circadian and activity prediction and probability of datapoint following the model;

Figure 15: Circadian model for a user;

Figure 16: Fitness as a function of q₀;

Figure 17: All factors risk R as a function of 0_O; and

Figure 18: Distribution of fitness and risk for the users in the database.

DESCRIPTION OF THE EMBODIMENTS

Before discussing the example embodiments in any more detail, first an overview will be provided.

A way of estimating resting heart rate is disclosed. This information may be used in the prediction of heart rate based on the moment of the day enabling one to determine if it is high or low. When estimating the resting heart rate a prior or base pattern for an average person is created from data measured from multiple users and a user’s heart rate is measured and fitted to this base pattern. The base pattern is formed of a cosinor with multiple, preferably four harmonics and a user’s measured heart rate is forced to fit this shape. The shape maybe horizontally translated to account for different sleep times of users such as night workers and students. A vertical translation can be applied to account for differences in resting heart rate due to fitness levels for example, and a scale factor indicative of an amplitude difference in heart rates between night and day may also be applied. Otherwise any user’s determined shape should match that of the base pattern. It should be noted that the vertical translation representative of the resting heart rate and the amplitude being representative of the difference between night and day both correlate with the fitness level of the user, that is a lower heart rate and larger range usually means better health status.

In summary expected HR in terms of circadian rhythm and physical activity is modelled. To account for differences between users latent variables are introduced that describe the general wellbeing of the user. The first part of the model predicts the user resting heart rate as a function of time (circadian rhythm), accounting for differences between users. The model will also take into account the change in heart rate caused by physical activity. The complete model is able to predict the user’s heart rate for any time of the day and activity level. Using this model actual HR values can be analysed and the probability for each HR data point to the model, be sensor noise, be model noise (caused by factors not part of the model) can be estimated.

Figure 1 shows a flow diagram illustrating steps in a method for generating the base pattern of an average user’s resting heart rate over a 24 hour period. In this method heart rate data is collected from multiple users over multiple days. The data from each user is analysed in unit time periods and it is determined if during each time period the user is estimated to be active or not. This information may be derivable from a wearable device configured to monitor steps for example or it maybe derived from changes in the heart rate data itself. If it is determined that the user is not active during this time period, the data is included in the sampled data. If it is estimated that the user is active, the data is discarded. Once the user’s data has been analysed for the whole 24 hour period, it is determined if the number of samples for the user is larger than a predetermined number and if it is the user’s data is added to the data for the generation of the prior. If it is not the user’s data is discarded. A subsequent user’s data is then analysed until all of the multiple users have been analysed. The collected data is then weighted so that lower heart rate data is given a higher weight than higher heart data as the lower heart rate data is more likely to be characteristic of a resting heart rate. The user’s data is then combined and fitted to a cosinor with multiple harmonics to generate a prior or base pattern.

This base pattern is then used in the analysis of individual’s data to estimate an individual’s resting heart rate as is shown schematically in figure 2 for example.

Figure 2 shows how the base pattern of Figure 1 might be used to estimate a user’s resting heart rate. Thus, a user’s heart rate is monitored from a wearable device. The user will enter their age and sex and the base pattern generated by the method of Figure 1 will be adjusted to account for differences that it is assessed this will make to the average user’s heart rate pattern. The user’s heart rate is then fitted to the adjusted base pattern and a model of the user’s heart rate is generated. From this, an indication of the average resting heart rate for the user can be determined and output.

Subsequent steps may then be performed where the resting heart rate for the user at a particular time is derived from the model and compared with the measured heart rate of the user and the difference is output or used to estimate a current work rate of the user for example.

One example of further steps that may be performed comprises estimating a likelihood for each measure, that is if the measured heart rate is close to the expected heart rate then the likelihood of that measure is high, whereas if it is remote from this value the likelihood is low. Once a set of likelihood values have been determined an average of these values during a time period is found, in some embodiments using the geometric mean. This average value may then be compared with an expected value for a corresponding time window and where it is determined to be different to the expected value by more than a predetermined amount an alert may be triggered. This alert may indicate abnormal physiological status, such as a stressful day, stressful meeting, flu, sleep disruptions or other things. In some cases an intervention may be triggered by the alert. The time window maybe 30, 40, 50... minutes long up to a time window of 6 hours. In some embodiments, time windows of increasing size may be used. So that an initially small window with low computational overhead may be used and where this indicates a potential problem, this may trigger further analysis in a longer time window. It has been found that the time averaged likelihood of a measured value is a particularly accurate indicator of an abnormal physiological status. Further details of a specific example of the generation of a base pattern or prior for resting heart rate are provided below.

Data is collected and gives as continuous HR and HRV (hear rate variable) data from almost 600 users. Activity data is collected for only a subset of those users.

Users that don't have at least 100 HR samples are discarded, so from the original set 224 users pass this filter, which is 91 percent of our user base. In summary we have a total of 247 users of which 160 are male (65% of the total) and 70 female.

Figures 3a and 3b provide a graphical representation of the different users in the sample group.

Circadian Rhythm

For each user the accelerometer activity was analysed and activities were labelled as Walking, Automotive, Still, and other labels. This is to separate data collected while moving, where HR is unlikely to reflect the true resting HR, and data collected while the user was still, or typing, or commuting. Activities were separated in 2 classes: Active (that includes the Walking label ), Inactive (that includes all other labels).

Prior for circadian rhythm

A single cosinor is generally used to assess the effect of circadian rhythm on heart rate. However, the shape of the HR as a function of the hour of the day significantly differs from a simple cosine: during the night the HR becomes slower for approximately 6 hours, and during the day it rises for approximately 18 hours, with 3 peaks

corresponding to meals. We are interested in describing the change in resting heart rate throughout the day and night, a model with a single cosinor would not allow us to properly estimate the expected resting heart rate, introducing a large error. We extend the cosinor model to use 4 harmonics, with a period of: 1 day; half a day; a quarter of a day; and an eighth of a day.

To build a prior for the resting heart rate we proceed as follows:

1. for each user we only keep the RR (resting heart rate) collected while Inactive

2. we sample 100 RR intervals, making sure we sample uniformly across the hour of the day

3. we only retain users that have at least 100 RR samples 4. we fit a periodic model for the test data using a cosinor with 4 harmonics. The dataset includes 74342 HR samples from 247 users.

The circadian model

We model HR over time with the sum of 4 sinusoidal components, as defined in the equation 4

where:

b is the average HR

the a parameters are the amplitudes of the harmonics

the f parameters are the phases of the harmonics

t is the time, expressed in hours

To fit the model we use harmonic regression.

We use data from all users, discarding individual differences, because we know that all our users have approximately the same habits (they work in the same office). This is a strong assumption and may introduce error, especially in high frequencies. The prior and posteriors are fitted at the same time, to account for users differences in phases. For the first iteration of the model, to keep the model simple, we train the prior before fitting the posterior.

Y(ί) = /3+a sin x+b cosx+csin 2x+c/ cos 2x+e sin 3x+ ^rcos 3x+g(sin 4x+/7 COs 4x

(2)

Because we expect some HR measure to be influenced by unknown factors such as mental stress, flue, that will result in a higher HR than the baseline, we use robust linear regression, setting weight equal to w(HR) = 1 . This weight HR function will give more importance to points with low HR.

The following relationships hold:

«i = a² + b²

b = 73.3942709707105

cr, = 4.6592966704931

• a₂ = 3.68637741057772

• cr₃ = 1.53797238362783

• cr₄ = 1.87123726362293

^ = -2.51687275875175

f₂ = 3.03496656460619

• <p₃ = 1.8336799004835

• f₄ = -0.103285661623034

Figure 4 shows the datapoints used for training and the resulting prior. The shape of the prior captures the shape commonly found in literature for the effect of circadian rhythm on HR.

Fit the model on the user

We use the above model as the prototype (and prior) for user’s HR throughout the day (circadian rhythm). However, every user has a different HR at night, during the day, and might go do sleep and wake up at different times than most of the users (that generated the model). For this reason we define a refined version of the model where some parameters are let free to adapt to user's data.

(3)

More concisely written as:

where: _• q₀ is the average HR of the user;

• q is the difference in HR range from the prior. Values larger than l indicate wider HR range, lower than l indicate narrower range;

q₂ is the offset in phase with respect to the prior.

This model forces the shape of the circadian rhythm to be the same for all users, the same as the prior. The only allowed changes are the average of the resting heart rate (0o), the amplitude of the oscillation between night and day (0i), and at what time the user goes to sleep and wakes up (the phase, 0₂).

We define a cost function J, shown in equation 5.

Where f_ί is the model of equation 3 refined with user parameters, applied to all the measures (m), and y is the actual measured value. Y(ΐ) is the function of the prior for the circadian baseline.

The first term is similar to R². However, because we expect external influences to make HR higher than expected (e.g. mental stress), we give more importance to y with low HR. This will make the model try to lower y instead of finding the average value.

The second term is a regularization term that adds a cost when f is far from the prior, behaving as a bayesian prior.

We have found an analytic solution to find the Q that minimizes J (see dedicated report at end of document). The computational cost of this operation is 0(m).

Now let's look at a particular user and see how the prior (black line) is adapted to the user data to a posterior model (red line) (see Figure 5).

Looking at a few more users in Figure 6 we see that the model generally adapts to the user data. In some cases the fit is not ideal. This will be addressed in future iterations of the algorithm.

Analysing all users Looking at the distribution of Q we notice that biases are present. This indicates that the prior is not entirely representative of the average user and being able to recognise this allows a future version of the algorithm to be improved.

Prior by age and sex

Figure 7 shows the distribution of Q by sex.

As expected from medical literature, average HR is slightly higher for females. Medical literature reports 5 as the expected difference, in our data the mean of Q for males is 72.62 and the mean of Q for females is 73.69, a difference of 1.07.

There are no differences on qi and Q2 by sex.

Correlation between age and qo is -0.08.

Correlation between age and Q 1 is -0.06.

This contradicts medical literature, which reports reduced circadian rhythm with increased age.

Modelling physical activity

Using the accelerometer on the band we estimated the activity of the user, with a granularity of 1 minute. We then define several metrics, counting the minutes active (active = not still) in the last x minutes. In particular we focus on minutes active in the last 5 minutes (to capture moments when the user is actively moving) and minutes active in the last 20 minutes (to capture sustained activity).

Prior for residuals by activity level

Figure 8 shows how many datapoints we have for each combination of activity values, in terms of number of minutes active in the last 5 minutes, and in the last 30 minutes. Because of the large disparity in number of datapoints, we show the log. We can see that for high values of activity in the last 20 minutes we have a small number of datapoints, statistics on that area should be handled with caution.

Figure 9 shows the average residuals (as relative error with respect to the actual HR) for each combination of activity levels. The average residuals increase gradually from o to 0.13 as activity increases. Table 1 shows the values.

Figure 10 shows the standard deviation of average residuals for each combination of activity levels. If we discard the outlier for o minutes in the last 5 minutes and 13 minutes in the last 20 minutes, the standard deviation gradually increases from 0.11 to 0.18 as activity levels increase. Table 2 shows the values.

Figure 11 shows the approximation of residuals using a linear model.

0 1 2 3 4 5

0 0.02

1 0.00 0 .02

2 0.00 0 .01 0.05

3 0.00 0.01 0.04 0.08

4 0.00 0.01 0.03 0.05 0.10

5 0.01 0.02 0.03 0.05 0.09 0.12

6 0.01 0.02 0.03 0.04 0.08 0.11

7 0.01 0.02 0.03 0.05 0.07 0.11

8 0.01 0.02 0.03 0.05 0.07 0.10

9 0.02 0.03 0.04 0.06 0.07 0.11

10 0.02 0.03 0.04 0.05 0.07 0.11

11 0.02 0.04 0.05 0.06 0.07 0.11

12 0.00 0.04 0.05 0.06 0.07 0.10

13 0.03 0.04 0.04 0.06 0.08 0.11

14 0.03 0.04 0.06 0.07 0.08 0.11

15 0.04 0.05 0.06 0.07 0.08 0.11

16 0.04 0.07 0.07 0.09 0.11

17 0.17 0.07 0.09 0.09 0.12

18 0.03 0.09 0.10 0.12

19 0.07 0.10 0.13

20 0.13 0.13

Table 1 : Relative residual means by minutes active in the last 5 and 20 minutes

User model with circadian and activity

We can now model the expected heart rate of the user as a stochastic function, de ned as the combination of circadian rhythm (t) and the physical activity (a). HR!J a) = (ί)( 1 + X(a)) (6) where X(b) ~ N(p(a), o(a)²) with m and a depending on the current user activity a, using the tables previously reported in this section. Figure 12 shows the distribution of the activity parameters using weighted linear fitting with weight equal to l/RHR. Adapt the activity model to user’s data

In this section we evaluate adapting the model of residuals by activity to each user.

Figure 13 shows the distribution of the HR relative prediction error training a linear model, using activity as independent variables, for each user and predicting the HR (using a 50/50 split), against the Null Hypothesis that predicts the HR residual using the tables by minutes active previously shown (that are the same for all users).

0 1 2 3 4 5

0 0.12

1 0.11 0.11

2 0.10 0.11 0.12

3 0.10 0.10 0.1 1 0.12

4 0.11 0.11 0.11 0.12 0.13

5 0.11 0.11 0.1 1 0.12 0.12 0.12

6 0.11 0.11 0.11 0.12 0.13 0.13

7 0.11 0.11 0.1 1 0.12 0.12 0.13

8 0.11 0.11 0.11 0.12 0.13 0.13

9 0.11 0.11 0.12 0.12 0.13 0.14

10 0.11 0.12 0.12 0.12 0.13 0.13

11 0.12 0.11 0.12 0.12 0.13 0.14

12 0.42 0.12 0.12 0.13 0.13 0.14

13 0.11 0.12 0.12 0.13 0.14 0.14

14 0.12 0.12 0.13 0.13 0.14 0.14

15 0.09 0.12 0.12 0.13 0.14 0.14

16 0.13 0.13 0.14 0.14 0.15

17 0.13 0.14 0.15 0.16

18 0.01 0.14 0.16 0.16

19 0.04 0.15 0.17

20 0.18

Table 2: Relative residual standard deviation by minutes active in the last 5 and 20 minutes

This experiment shows that there is no significant increase in precision in adapting the prediction model on every user (personalized prediction model relative error mean 0.0157278989303904, with standard deviation 0.121327599572675, against Null

Hypothesis relative error mean 0.0041268686240701 and standard deviation

0.124722773692561). Therefore we will not adapt the activity model around each user, and we'll use the residual prediction by activity level tables previously shown. However, fitting a model for each user is interesting to analyse the response of the user to physical activity, to estimate wellbeing and fitness. Probability of outliers

We can now, for each HR measure, calculate the probability that the datapoint is in line with the model we created.

Figure 14 shows 2 graphs. The top graph shows the actual measured HR (black dots), circadian baseline (green line), predicted HR (red dots) with un-certainty (error bars). The bottom graph shows, for every datapoint in the top graph, the probability that it is generated by the model.

It is now possible to identify outliers (single datapoints with low probability) and sections not explained by the model, therefore interesting because they show abnormal HR that could be caused by other factors such as mental stress.

General wellbeing

In several studies the correlation between resting heart rate (RHR) and physical fitness is analysed, and its correlation with increased risk of mortality.

As shown in Figure 15, RHR can be estimated from qq, that expresses the average RHR throughout the 24 hours. The Maximum RHR is approximately equal to qq + 5Q1 and the minimum RHR is approximately qq Iqqi. During the day RHR oscillates between qo and 0o + 5 0i.

We can de nearest approximation of the fitness of the user using 0o and the conclusions from previous studies.

The fitness function F is defined as a sigmoid function with the flex at 70, as shown in equation 7 and Figure 16

As previously discussed, the fitness function F is designed to report high values for 0_O below 50, normal values for 0_O between 50 and 90, and low values for 0_O above 90.

Previous studies have concluded that“risk of mortality increases with 16% per 10 beats per minute (bpm)". RHR less than 5obpm is correlated with maximum physical fitness, RHR more than 9obpm with minimum fitness. As shown in equation 8 and Figure 17, the mortality risk function R is defined as a function of Q : max(5o.min(9o.9₀))-5o

R = 30 . 1.16 (8)

Figure 18 shows the distributions of the estimated user fitness, using the function F , and the estimated all-factors mortality risk, using the function R.

A further advantage of embodiments is that even with data that contains missing datapoints and is unevenly sampled, the mathematical tools are such as to enable the data that is not evenly sampled to be fitted to the model. This is implemented at a mathematical level which is described in the following Periodic Pattern Fitting summary. It's an approach similar to a Fourier Transform (it works in the domain of frequency instead of time), but without the restrictions and assumptions of traditional Fourier analysis. Periodic Pattern Fitting describes the problem that it solves starting from an abstract/high level approach, section 1, and considers special cases, every time more specialized, down to our special case (we want to fit our datapoints to a set of periodic functions, of known period, only allowing a limited set of parameters to change, thus preserving the "shape"). Section 2 explains how this is implemented numerically, and why it's computationally efficient.

1 Problem definition

let /*(*) t C 4 S* is a. fa ly of toctlone a&d (_.¾¾^: _5%¾}o<_{ί' '}>_{> &} sol· of weighted data pohto (a? f C, f f S, f g ^T ) We waat to tied the best parameters « ~ {«}_j? ^v {*?}

tha

as close as pfsemble to the set of data po te, that fe it imfees ties square loss

This k as always postive real fuette* of the t ¾cw « \4 and thus to mrnxmam has its total idereefiu! o< ¾l t ¾m that is

whets a¹ deato# the

Giw that ¾ d w>t do ant oa tto moatioua! pu oeteis mtd v to iofai differential is ¾ao_: thus WO 000 p ih e

*8 a ector space o*r C where a* enotes the co plex conjugation of a

We thus rtssifnte the expanhon of the t o total differentials 4_nxV(¾-) and

The input data appears only hr tie; hmer nx m (over )_t we cars ti e rewrite these, «sing the .«otmtk ::

rb fLnf iiii

The e uation sysieta in Eqa. 3( becomes then (* hns been replaced by k for readability}:

(40)

Which is a linear syste «5: both P au·,! S are not dependent on v.

1,2 Case 2; periodic fooetlons

A more general case is w hen the base fun·'^; .au*. fm the fit are generic eriodic fimetlons whose period b nndf iple of fire sann_* b; u e period, La.:

o (bat

then Eq. 26,21,28,20 her'e es

using «¾ ··- e^iH¾',¾ _t and

similarly for _«F:

1,2,1 Si ple sha e 81

parifonlar of the ab ve is wkmi there «re only 2 base perfodk fuiioicti* a coostaBt function /o(as) 1, and an arbitrary m awragp periodic fooeikm fi{x) Liy ^< e Fί *^'^ _' In Ibis ease there are only three free para eters:

* ¾% «% which is a vertical tranelalion of the fitting: shape /_{¾ ···} /

* ¾¾ - fo hich la a scale factor for /

* w; ¾ which is a» hor oiifat tanslation of the lif ing shape.

The a px i afmg function is fhns:

the parameters are thus:

Also, since there is no U_{Q '}variable we have

BUQ · · u {61

We sate that the camiaon factor m equations of 55 sim lifies as;

- fyfy, - å j^{. (}¼, .._·/ Jfy (65) thus we can rewrite Eqp 56 (assuming, to simplify somm f — 0} as

Were ioofcmg for solmk s where 4 # Cl beoroae otherwise the tel pe” we are

fitting 5a Hot· presmtt_* also _s% L The equation system becomes

fo thte reason, tdtipiy g both aides of the sx acfo» by z^^~i ( Meh e know is not aero, sir»» <) m an exjxaientSai) wo obtain a ohme ial of degree 6c - 2. 2 Numerical solution

A way to solve the polynomial

method. The com anion matrix of a polynomial is a sparse matrix which can he crafted to have the etiar&eterisfte polynomial which is 0 (moMe) multiple of the original olvoorfiia!². In this way, the mafrixfo eigenvalues are also the root of the original polynomial. This equivalence allow finding approximately the roots of the polynomial using: any numerical analysis algorithm on the matrix (see https i //sn . «ikipsdla, oig/aiki/ ist_0f i »sricai_anaiysis to_jpics#Iign»valna algorithms for a list)

2,1 Complexity and engineering considerations

This method of fin ing the best fit for a periodic shape has some interesting points to consider;

1 the method find the best approximation of data distributed In the x axis

without any regularity, with any weight istribution

2, the data dependency is only present in the final equations through the

“s ape coefficients⁵' S4, which depend only on ay, ami the“limner convolution coefficients" which depend on both a¾ ami Both these co fick ts cm be computed independently for every single point of data, and then aggregated together by averaging them. This allow an kficiont distribution of the computation, which cm even be performed online f n?., by the device themselves as data are aerpired) Moreover, ehaugmg shape do not alftet compntatkni of these ciMilfkiente , winch Is the most com putationally expensive step In the analysis (see below):: this allows fitliug

any shapes with the saute data with minimal overhead (see below), which is useful in peak detection

3 both methods shown (low-component trigonometric lit and periodic shape

tit) have lo complexity: compiling the tables of f¾ and i¾ is linear in the numbe of data points an in the reeolntion teu hef of eompom.nda required to describe the shape). The following steps which are signitlcant

In term of complexity are:

* In the ease of trigonometric fit solving a linear system where the

number of equations 1» the number of trigmm eifls fmtethm used

* In flic case of shape ft finding the eigenvalue of a sqn c matrix

whose m is the number of foarkr component time describe the .dra

but elss giiig the * period arid the nu ber of feircmwt s <lo, «» to sns wh ii a both i ependent on the iimhber of data points».

3 Appendix: c^A^? _V(x) = 4_UV(#)

Since

we nan write the total differential of it by splkafcty separating the real and imaginary part of the variable, that Is:

4 Unsorted

temp: t is shcmid noi hv fteg i -n

Abstract— Wrist-worn wearable devices equipped with heart rate sensors have become increasingly popular. The ability to correctly interpret the collected data is fundamental to analyse users well-being and perform early detection of abnormal physiological data. Circadian rhythm is a strong factor of variability in heart rate, yet few models attempt to accurately model its effect on heart rate.

We present a novel mathematical model that allows the accuracy of multiple component cosinor model. At the same time the parameters of the model are easy to interpret parameters: MESOR, Amplitude and Acrophase, the same produced by the single component cosinor model. Moreover, our model does not require uniformly sampled data, and allows to weight each datapoint independently.

The model is computationally cheap, depending linearly on the size of the data. The computation of the model does not need the full dataset, but only two surrogates ( F and S ), where the data is accumulated. This implies that the model can be implemented in a streaming approach, where the data is accumulated on the device, with important consequences for security and privacy of the data, which never leaves the user devices. We show the accuracy of our model comparing its performances against the single component cosinor model by measuring HR prediction RMSE. We also show that the model parameters can be used to detect sleep disruption in a qualitative experiment.

I. INTRODUCTION

In the last decade, wearable activity trackers have become more and more popular to monitor Heart Rate (HR). Being able to predict a user resting HR during the day is useful to interpret data collected from wearable devices and to make inference about people health status. As a matter of fact the resting HR is widely investigated in literature because it is an independent predictor of cardiovascular and all-cause mortality in men and women.

Physiological measures such as HR are continuous variables that are subject to the influence of external and internal stimulus changing over a daily cycle of about 24h with patterns that have been defined as circadian rhythms. The time series obtained from wearable device are often noisy (i.e., variations in the biological system that are not part of the deterministic portion of the signal and that are derived from external errors such as instrument inaccuracy), short (i.e., low sampling) and sparse (i.e., unequal time intervals between observations). Despite these time series problems and the fact that the real signal is not a simple sinusoidal, HR circadian rhythms analysis is usually done using the cosinor method with a single component because the resulting parameters can be easily interpreted. In particular, the cosinor curve -with single component provides 3 parameters that describe the circadian rhythm: i) Mid-line Estimating Statistic Of Rhythm (MESOR) is an estimation of central tendency of the distribution of values across the cycles of the circadian rhythm computed using a cosine function; ii) Amplitude is the difference between the peak and the mean value of a wave; iii) Acrophase is the time of the day at which the peak of the circadian rhythm occurs.

It was found that the shape of the HR as a function of the hour of the day significantly differs from a simple cosine: during the night the HR becomes slower for approximately 6 hours, while during the day it rises for approximately 18 hours, with 3 peaks corresponding to meals [5], [6]. However, only a few studies in literature propose a cosinor method with multiple-component to analyze periodic components in non- sinusoidal longitudinal time series [3], [4]. The results of the multiple-component cosinor analysis include several parameters which are not easy to be interpreted in comparison with the three parameters provided by single-component cosinor analysis. Despite the multiple-component model was found to be more suitable to approximate the signals waveform when it deviates from sinusoidal, the single-component cosinor model is still widely used to investigate physiological circadian rhythm because it provides easy-to-interpret parameters [4].

To the best of our knowledge, this is the first study that provides an algorithm able to predict the resting HR during (RHR) the day correcting the actual HR by the effect of external stimulus such as physical activity and stress status. In particular, we are interested in describing the change in resting HR throughout the day and night by using a multiple component cosinor where a model with a single model would not allow us to properly estimate the expected resting HR because it would introduce a large error. Hence, the aim of this study is to define a non-sinusoidal model in order to assess the user’s resting HR on the base of the population resting HR circadian rhythms, while providing parameters that are easy to interpret.

A. Contributions of This Paper

In this paper we present a mathematical derivation of the single-component cosinor model with multiple components that fits an arbitrary function (i.e., the model which describes the resting HR circadian rhythm of the population) on a user data, thus permitting us to predict the user’s resting HR with a high accuracy. Fitting the user’s data on a arbitrary function makes it possible to use weighted data derived from the population distribution to detect resting HR by using noisy time series (e.g., data derived from unreliable source such as wearable devices) without worrying about outliers and resulting skewed models.

Unlike the cosinor model with multiple components provided by Cornelissen, the model proposed in our study provides the same three easy-to-interpret parameters which can be derived from the cosinor model with single component (i.e., MESOR, Amplitude, Acrophase ).

Our model allows the assignment of weights to every data point used for fitting the expected shape of to the user data. This allows the implementation of advanced analyses like assigning less importance to data points collected while the user was moving (expecting motion artefacts). Using weights also makes it possible to implement real-time streaming analysis of the data, performing a new analysis as new data arrives, assigning less importance to older data.

The model presented in this paper is computationally light and, for this reason, it is suitable to be executed on a wearable device. The complexity of the accumulation phase is linear on the number of data points, while the complexity resolution phase depends on the number of frequencies used and not from the number of data points. II. BACKGROUND

Short and sparse time series were historically analyzed by the single-component cosinor model. However, this model is not suitable for noisy time series (i.e., the one showing several outlier values) such as the one recorded by wearable activity trackers. Cornelissen at al. [4] easily extends the single-component cosinor to a multiple- component model in order to analyze the time series in chronobiology. Instead of solving a system of 3 equations in 3 unknowns to define the parameters of the circadian rhythm curve (i.e., MESOR, Amplitude, Acrophase ) there are 2c + 1 normal equations to estimate MESOR and p pairs of Amplitude and Acrophase (i.e., c refers to the number of harmonics selected to describe the circadian). This approach better approximates the signals waveform of the circadian rhythm compared to single component cosinor model thanks to the higher number of cosinor components used to fit the model on the circadian rhythm curve. For example, a 2 component cosinor model that describes periods of 24 and 12 hours has been found to be well-suited to approximate the nightly drop of blood pressure [8]. However, the parameters provided by this model result hard to interpret and to compare with the results provided from single- or other multiple component cosinor methods [3], e.g. a small change in the phase of an harmonic will result in drastic changes in the resulting shape of the function.

III. METHODS

A. The mathematical model

1) Prior from general population: We model the prior of RHR over time A K) with the sum of c cosines, as defined in the multiple-component cosinor in Cornelissen at al. [4]. This equation can be rewritten using Euler’s complex exponential formula as Equation 1. The number c of components to use is a free parameter.

The parameters b, a, and 0 correspond to the MESOR, Amplitude, and Acrophase, and can be found as explained in Cornelissen at al. [4]. If we define y₀ =

can rewrite Equation l as Equation 2

k=-c

Note that we have defined y so that ip__k is equal to the complex conjugate of xp_k . Therefore, the number of real parameters has not changed.

2) Fitting user data with locked parameters: Single component cosinor can only accurately fit data modulated by a simple sinusoidal factor. Circadian effect on RHR is unlikely to be fully captured by such a simple function, therefore we expect that using more components we will be able to fit the data more accurately. Multiple component cosinor can fit data more accurately [4], but the resulting set of fitting parameters, (a MESOR, Acrophase, and Amplitude for each component) can not be easily interpreted, i.e. a small change in the phase of one of the components dramatically changes the resulting shape, making it difficult to interpret the parameters of each component. Moreover, fitting a model for each user with all the free parameters will fit shapes that might not be what circadian is supposed to be, and it will be difficult to get an idea of what is the phase and amplitude.

We use the A(h) model provided in the previous section as the prototype (and prior) for users RHR throughout the day (circadian rhythm). However, every user has a different RHR at night, during the day, and might go do sleep and wake up at different times than most of the users (that generated the model). For this reason we define a refined version of the model where some parameters are let free to adapt to users data. In particular, Equation 3 shows the formula that locks parameters of the c components only allowing 3 degrees of freedom (the same used in single component cosinor analysis). This is equivalent to single component cosinor analysis, using a different function than a cosine.

The free parameters in Equation 3 are:

q₀ the average RHR of the user;

0_! the amplitude of the oscillation of the user RHR with respect to the prior.

Values of 0_! close to 1 indicate that the range of the resting heart rate is similar to the prior, less than 1 indicate that the user resting heart rate has a smaller range, etc.;

q₂ the phase of the user circadian rhythm with respect to the phase of a cosine with period equal to one day. Values of q₂ close to o indicate that the user is sleeping at midday. q₂ ranges between -0.5 and 0.5. The Q parameters allow to express the user RHR, U (h), as a vertical shift ( q₀ ), horizontal shift (q₂ ), and a change in amplitude (0_! ) of the population prior A(h).

3) Loss function: To find Q that best approximate U(h) for a particular user, we fit U (h) to the RHR data collected for that user, minimising the distance L between the data and U(h) as shown in Equation 4, where y₇- is the y^th measure of RHR of the user, x_j is the time of the measure, and w₇- is the weight assigned to the ;^th measure.

Finding the Q that minimises L is not trivial, because we use periodic functions that are not simple cosines, but a particular combination of cosines, defined by the cx and 0 parameters, that are not allowed to change. The literature of cosinor based rhythmometiy does not offer any solutions, neither analytic nor numeric. Therefore, we developed an analytic solution to finding the Q that minimises L.

First, we introduce the following notation to simplify further formulas:

We can now develop Equation 4 using the definition of U in Equation 3 into the following:

i(0) = ^ w_j yj + 0o ^ w_} + Q² ^ 0_m 0_k <;^k+mS_k+m

] j k,m

= -20₀ Fo + 0₁ å 0m C,^m (20₀ S_m -2 F_m ) (8)

m

where the indices k and m range between— c and c.

To minimise the loss we put to zero the partial derivatives with respect to q₀ , q₁ , q₂ we obtain the system in Equation 9

Substituting 0_O and then Q₁ in the equations of the system 9 we find a polynomial in :

Since the indices m, L and k range between— c and c, and Vm.

= 0, the polynomial in Equation 10 has degree 6c -2.

The roots of the polynomial in can be found with spectral methods, i.e. using the eigenvalues of the companion matrix [9], and then refining the solutions from numerical error using the Newton- Raphson method.

From the roots of the polynomial we can determinate the values of q₂, and subsequently of q_± and q₀, where the total derivative of 1(6) is zero. We compute (6) using Equation 8 on these points to determine the absolute minimum.

4) Engineering considerations: The system in Equation 9 depends on data indirectly through the values of F and S. Moreover, also the square loss function can be expressed in term of these values (and the total sum of yf). This means the whole algorithm does not require access to the full dataset, but only a limited set of values (namely the values of S, F and the sum of squares of y values if you want to output an absolute value for the loss). All these values are moreover an additive function of data sets, that is for two disjoint set of values A = {x_t ,y_t }and B = {x_j , y_j } S _AUB = S _A + S_{k B} and F _AUB = F_{k A} + F_R,B · This implies that it can be computed separately for each subset of samples, and then added together to obtain the total value: for example it can be accumulated directly in the user’s device, with big advantages in computational complexity, data transfer and privacy. In particular the asymptotic analysis of the algorithm complexity is the following (where N is the number of data points, c the number of components):

• The accumulation of data into the surrogate values S and F is 0(N c)

• To solve the polynomial we need to find eigenvalues of a 6c -2 sided matrix which requires 0(c³) operations

• To compute the loss for a given set of fit parameters we required 0(c²) operations, which has to be repeated for each candidate solution, for a total of

0(c³) operations

The total complexity is thus 0(c³ + Nc), in particular it only depends linearly on the size of data.

B. Experimental setting

To test the accuracy of our model of circadian effect on RHR against the traditionally used single component consior, we used the data collected by Biobeats through their corporate wellbeing product Biobase, a stress management app that uses a wrist worn wearable device to passively collect activity data (steps), sleep phases, and heart rate. To test our model we used a dataset collected by Biobeats from two pilots with financial sector organizations. 47 individuals (F=52%, M=47%), with an age range of 20 - 63 years (M = 29.29, SD = 11.38), wore a wrist worn wearable device that collected heart rate every 10 minutes, activity data (steps) every 20 seconds, and sleep phases every night, for four weeks. The dataset has a total of 163179 HR measures.

1) Building a prior by fitting b and y on the general population: To build a prior for the resting HR we proceed as follows:

• for each user we only keep time series with at least 100 heart rate samples;

• we split the dataset using a leave on out strategy, training the prior on a dataset that did not include the user for which we wanted to calculate the Q parameters;

• for each HR measure, we calculate the amount of steps taken by the user in the previous 5 minutes;

• we remove from the dataset the HR points where the number of steps in the previous 5 minutes was larger than 10. This step is carried out to ensure we estimate RHR, instead of HR, that is influenced by the physiological request of physical activity;

• we perform cosinor analysis on a single component to find the Acrophase of this user (every user will generally have different sleep habits);

• we align the data of all users, forcing all users to have the same phase as a cosine with period equal to 1 day;

• we use multiple-components cosinor using c = 4 cosines, with period of 1 day, half a day, a quarter of a day and an eight of a day.

Figure 1 shows one of the priors fitted using the leave out-out approach. In that case the values for b and y, from equation 2, are:

• b = 68.75;

• fΐ = 1.63;

• y2 = 0.69e ^{i3 04};

· f3 = 0.13e ^{i0 19};

• f4 = 0.19e^_i1·⁵⁴;

b represents the population average. _± describes the first of the 4 components cosinor fitting. The modulus expresses how important is this component in the fitting, and the argument expresses the phase of this component. has argument equal to zero, because we explicitly aligned all users to a cosine. We can see that the first component captures 62% of the overall information, as the modulus of 1/q is equal to 0.62 multiplied by the sum of the moduli of all _±, f₂, f₃, and f₄.

2) Building the models of a user resting heart rate by fitting Q on the data of that user: We then selected the data of the target user. For every HR measure we calculated the surrogates S and F, solved the polynomial in V, finding the Q parameters of Equation 3 that minimised the loss function. We applied exponential decay to the previous data. The decay function is applied to the weights wi used in the loss function L. w_t = e^{~ t}, where At is the time, in days, between the datapoint i and the latest datapoint, that triggered the update of the model.

Figure 2 shows the data of one day of a user, with the fitted single component cosinor model and our circadian model. We can see that our model better captures the sudden change in HR when the user awakens, and the fact that the time awake is longer than the time asleep.

Figure 3 shows seven days of HR and sleep data of one of the users, and the estimated circadian parameters:

• the first graph shows the HR measures of the user as black dots, and the RHR as estimated by the circadian model as red dots;

• the second graph shows the estimated MESOR, q₀, evolving over time, changing for every HR measure;

• the third graph shows the estimated Amplitude, q₁, evolving over time, changing for every HR measure;

• the fourth graph shows the estimated Acrophase, q₂, evolving over time, changing for every HR measure;

• the fifth graph shows the sleep data measures by the wearable, in the graph we show at what time the user went to bed and woke up, as a line.

It can be seen that the estimation of the RHR of the user fits the user data better than a single cosine could do. Looking at the estimated circadian parameters, it can be seen that at the fourth night, slightly after timestamp=i546300000, the MESOR increases and the Amplitude decreases, indicating non optimal physical conditions of the user, and that the Acrophase suddenly changes, increasing from the stable average. Observing the sleep data from the band it can be seen that on the fourth night the user suddenly changed sleep habits (from a healthy stable time to bed of about 11pm, to 3am). The lack of proper sleep can be seen in the HR data and is reflected in the changes in 9₀ and q_±. The sudden change in the user sleep habit is reflected in the sudden change in q₂·

3) Predict user’s resting heart rate: We applied the method described in the previous sections to obtain an evolving model of the users. We used those models to predict the next HR, given all the user’s previous HR datapoints. We compared the prediction accuracy of our model against the single component cosinor model.

The cosinor RMSE is 5.73 ± 2.33, our circadian model RMSE is 5.15 ± 1.04. Our circadian model RMSE is in average 10% lower than the cosinor RMSE. A paired t-test returns, with a 95% confidence, a difference between cosinor RMSE and our circadian RMSE between 0.04 and 1.13 (p-value< 0.05).

We also tested our model against the Null Hypothesis that simply using the previous HR measure to predict the current HR measure yields better results than modelling circadian effect on HR. The Null Hypothesis model has RMSE equal to 5.88±I.86.

The loss function L{9) can be used to calculate the standard error, useful to have an estimation of the distribution of the predicted data. This information could be used to automatically analyse HR activity, finding anomalous data points that could be caused by factors not modelled by circadian effect on RHR, such as physiological conditions.

IV. CONCLUSION

In this paper we presented a novel mathematical model that can be used to approximate a user’s resting heart rate. The model achieves the expressiveness of the multiple component cosinor, i.e. the model fits the data. At the same time the parameters of the model are easy to interpret parameters: MESOR, Amplitude and Acrophase, the same produced by the single component cosinor model. The model does not require uniformly sampled data, and allows each datapoint to be weighted independently.

The model is computationally cheap, depending linearly on the size of the data. The computation of the model does not need the full dataset, but only two surrogates (F and S), where the data is accumulated. This implies that the model can be implemented in a streaming approach, where the data is accumulated on the device, with important consequences for security and privacy of the data, which never leaves the user devices. The model is, to the best of our knowledge, the only computational model that can be used to predict users resting heart rate and to analyse the data using the parameters traditionally used to describe the circadian modulation of physiological activity.

The model we developed goes beyond fitting circadian activity on resting heart rate, and it can be used to fit arbitrary periodic real valued time series, but also vectorial data (e.g. gesture recognition from the accelerometer), or complex data. With an extension to the mathematical model it could be used to fit non periodic data using a linear combination of an arbitrary set of functions.

A person of skill in the art would readily recognize that steps of various above- described methods can be performed by programmed computers. Herein, some embodiments are also intended to cover program storage devices, e.g., digital data storage media, which are machine or computer readable and encode machine- executable or computer-executable programs of instructions, wherein said instructions perform some or all of the steps of said above-described methods. The program storage devices maybe, e.g., digital memories, magnetic storage media such as a magnetic disks and magnetic tapes, hard drives, or optically readable digital data storage media. The embodiments are also intended to cover computers programmed to perform said steps of the above-described methods.

Although embodiments of the present invention have been described in the preceding paragraphs with reference to various examples, it should be appreciated that modifications to the examples given can be made without departing from the scope of the invention as claimed

Features described in the preceding description may be used in combinations other than the combinations explicitly described.

Although illustrative embodiments of the invention have been disclosed in detail herein, with reference to the accompanying drawings, it is understood that the invention is not limited to the precise embodiment and that various changes and modifications can be effected therein by one skilled in the art without departing from the scope of the invention as defined by the appended claims and their equivalents.

Claims

1. A method of determining a user’s resting heart rate, comprising:

receiving data indicative of a user’s heart rate sensed over a period of time; fitting at least some of said sensed data to a base pattern indicative of estimated variations in an average person’s resting heart rate over a predetermined time period, said base pattern comprising a cosinor comprising a plurality of harmonics; and

determining said resting heart rate from said sensed data fitted to said base pattern.

2. A method according to claim l, wherein said predetermined time period comprises 24 hours.

3. A method according to claim 2, wherein said plurality of harmonics comprise four harmonics, said four harmonics comprising one day, a half day, a quarter of a day and an eighth of a day.

4. A method according to any preceding claim, wherein said sampling is performed over several days and said method further comprises the step of combining heart rate data for a same period of time from different days in said sampling period.

5. A method according to any preceding claim, further comprising determining at least one of sex and age of said user and adjusting said base pattern in dependence upon said at least one of said sex and age of said user prior to performing said fitting step.

6. A method according to any preceding claim, further comprising determining said user’s activity levels and selecting heart rate data during periods of time where activity level is low and discarding heart rate data during periods of time where activity level is deemed to be high.

7. A method according to any preceding claim, comprising applying a weighting to at least some of said user data prior to performing said combining step.

8. A method according to claim 7, wherein said step of applying said weighting to said data comprises applying a lower weighting to higher heart rate data than to lower heart rate data.

9. A method according to claim 7 or 8, wherein said step of applying said weighting comprises applying a higher weighting to newer more recently measured data than to older data.

10. A method according to any preceding claim, wherein said step of fitting comprises forcing said user’s data to have a same basic shape as said base pattern, said same basic shape having a freedom of a vertical translation, a scale factor and a horizontal translation with respect to said base pattern.

11. A method according to any preceding claim, wherein said step of fitting comprises forcing said user’s data to have a same basic shape as said base pattern by performing harmonic regression on said data to force the shape of the variation in heart rate overtime to be substantially the same as the base pattern, any variations between the base pattern and user’s pattern being in the scale factor of the shape the vertical translation of the shape , and in the horizontal translation or phase offset of the shape with respect to the base pattern.

12. A method according to claim 10 or 11, wherein said step of fitting comprises applying said at least some of said user data to the cosinor of the base pattern and determining a cost function by squaring the difference between measured heart rate and modelled heat rate weighted to favour lower heart rates, and determining an analytic solution that minimises the cost function.

13. A method according to claim 12, wherein said cost function to be minimised comprises a further regularisation term that adds a cost where the cosinor with the user’s data is far from the base pattern.

14. A method according to any one of claims 10 to 13, wherein the vertical translation of the basic shape to the base pattern, the scale factor and the horizontal translation provide an indication of average resting heart rate, amplitude difference in heart rate between night and day, and an indication of sleep times of the user.

15. A method according to any preceding claim, wherein said step of determining said resting heart comprises determining at least one of an average resting heart rate and a resting heart rate as a function of time of day.

16. A method according to any preceding claim, comprising a further step of outputting an indication of said resting heart rate, wherein said indication comprises at least one of an average value of a determined resting heart rate, a current resting heart rate and an indication of fitness.

17. A method according to any preceding claim, wherein said at least some of said data comprises two surrogates derived from a measure of the resting heart rate and the time of the measure.

18. A method according to claim 17, wherein said two surrogates comprise a shape coefficient S and Fourier convolution coefficient F.

19. A method according to any preceding claim, wherein said at least some of said data is stored on said device.

20. A method of comparing a user’s heart rate to a predicted value, said method comprising:

sensing a user’s activity level and determining an expected heart rate from an estimated increase in heart rate due to said sensed activity level and a determined resting heart rate for that time of day, said resting heart rate for said time of day being determined by a method according to any one of claims 1 to 18;

outputting an indication of a difference in measured heart rate and predicted heart rate.

21. A method according to claim 20, comprising determining from said difference in measured heart rate and predicted heart rate a likelihood of said measured heart rate occurring.

22 A method according to claim 21, comprising determining a plurality of said likelihoods of a plurality of measured heart rates occurring, said plurality of measured heart rates being measured over a predetermined time period, and determining an average of said plurality of determined likelihoods during said time period.

23. A method according to claim 22, wherein said average comprises a geometric mean of said likelihood values for said time period.

24. A method according to any one of claims 22 or 23, wherein where said averaged difference is greater than a predetermined value, said indication comprises an alert.

25. A method according to any one of claims 20 to 24, wherein said indication comprises a wellbeing indication.

26. A method of determining a base pattern indicative of resting heart rate variations in an average person over a period of time, comprising:

receiving data from a plurality of users indicating heart rate over said period of time;

selecting heart rate data during periods of time where activity level is deemed to be low and discarding heart rate data during periods of time where activity level is deemed to be high;

combining said data from multiple users and multiple days;

generating a base pattern indicative of estimated changes in resting heart rate over time for an average person by fitting a periodic model to said combined data, said periodic model comprising a cosinor with a plurality of harmonics.

27. A method according to claim 26, wherein said step of combining comprises an initial step of averaging heart rate data from each user for a same period of time in a day from a sample period of several days and coalescing said data for a plurality of users.

28. A method according to claim 26 or 27, wherein said step of coalescing comprises combining said data, said data being weighted prior to said combination, higher heart rate data having a lower weighting.

29. A method according to any one of claims 26 to 28, wherein said step of fitting comprises using harmonic regression to said sampled data to fit said sampled data to said shape defined by said cosinor.

30. A method according to any one of claims 26 to 29, wherein said plurality of harmonics comprise four harmonics, said four harmonics comprising one day, a half day, a quarter of a day and an eighth of a day.

31. A method of determining a user’s resting heart rate according to any one of claims 1- 25, wherein said base pattern is determined according to a method according to any one of claims 26 to 30.

32. A computer program comprising computer executable instructions which when executed by a processor are operable to control said processor to perform a method according to any one of claims 1 to 31.

33. A device for determining a user’s resting heart rate, comprising:

a sensor configured to continually sense said user’s heart rate over a period of time;

processing circuitry configured to fit at least some of said sensed data to a base pattern indicative of estimated variations in an average person’s resting heart rate over a predetermined time period, said base pattern being modelled by a cosinor comprising a plurality of harmonics and to determine said resting heart rate from said sensed data fitted to said base pattern.