### Study design and participants

Our study uses workout measurements contributed to the Apple Heart and Movement Study (AHMS) between November 2019 and July 2022^{18}. The Apple Heart and Movement Study is a prospective, single-group, open-label, siteless, pragmatic observational study conducted in collaboration with the American Heart Association and Brigham and Women’s Hospital to investigate the relationship between physical activity, mobility, and heart health. The study was approved by the Advarra Central Institutional Review Board (PRO00036784) and registered to ClinicalTrials.gov (ClinicalTrials.gov Identifier: NCT04198194). There is no compensation for participation.

The AHMS allows participants to securely share information collected by their Apple Watch. To be eligible, the participants must own an Apple Watch paired with their iPhone, be at least 18 years old (at least 19 years old in Alabama and Nebraska; at least 21 years old in Puerto Rico), live in the United States of America, have installed the Apple Research app on their iPhone, do not share their iCloud account, iPhone, or Apple Watch with anyone else, and are willing and able to provide informed consent to participate in the study. The study app was used to verify eligibility, obtain participants’ consent, provide study education, and direct participants through the study procedures.

Within the AHMS, we selected participants who chose to use the Workout app on their Apple Watch to record their outdoor runs. The AHMS logs the outdoor weather information *W* (temperature and humidity) at the time of each workout and provides the participant’s heart rate during each workout along with four measurements of exercise intensity: the instantaneous speed from the pedometer sensor, the instantaneous speed from the GPS, the step cadence from the pedometer, and the elevation gain from the altimeter. The sensor measurements are interpolated on a 10-s grid to form, for each workout *w*, a heart rate time-series \(\widehat{{{{\rm{HR}}}}}\in {{\mathbb{R}}}^{d}\) and a multivariate time-series of exercise intensity \(I\in {{\mathbb{R}}}^{4\times d}\); where *d* is the duration of the workout. We finally only selected workouts with a duration between 15 and 120 min long that had no missing data and further filtered out participants with less than 10 valid runs. We obtained 270,707 outdoor runs from 7465 distinct subjects (see Fig. 1 for cohort summary and inclusion criteria).

Our hybrid approach blends a physiological model of heart rate dynamics with machine learning components (e.g., deep neural networks) to adapt it to wearable data and personalize it to individual subjects. We first detail the physiological model and then describe our ML-based augmentations and learning algorithm. A diagram of the end-to-end system is detailed in Fig. 2.

### Physiological HR dynamics model

Several works in the sports physiology literature have studied heart rate dynamics in response to exercise using ordinary differential equations (ODEs)^{14,15,16,17}. These approaches translate the physical mechanisms of the human body into differential equations in order to incorporate domain (physiology) knowledge in the modeling. This is an appealing method to build interpretable and ultimately trustworthy models of fitness and health.

A common approach for modeling changes in heart rate, HR, due to exercise intensity, *t* ↦ *I*(*t*), is to introduce the oxygen demand, *D*, as an intermediary quantity through coupled ODEs^{15}

$$\left\{\begin{array}{ll}\dot{D}(t)\;\;\;\,=B\cdot \left(f(I(t))-D(t)\right),\\ \dot{{{{\rm{HR}}}}}(t)\;\,=A\cdot {\left({{{\rm{HR}}}}(t)-{{{{\rm{HR}}}}}_{{\mathrm{min}}}\right)}^{\alpha }\cdot {\left({{{{\rm{HR}}}}}_{{\mathrm{max}}}-{{{\rm{HR}}}}(t)\right)}^{\beta }\cdot \left(D(t)-{{{\rm{HR}}}}(t)\right),\\ {{{\rm{HR}}}}(0)\;={{{{\rm{HR}}}}}_{0},\\ D(0)\;\;\;={D}_{0}.\end{array}\right.$$

(1)

In this dynamical system, the function *f* translates the instantaneous activity intensity *I* into the necessary oxygen demand for *I*. The top equation attempts to match the current body oxygen demand *D* with the instantaneous demand *f*(*I*) (also known as the *drive function*). Parameter *B* controls how fast *D* adapts to *f*(*I*). At the same time, the second equation drives the heart rate HR toward the pace required to deliver the demand *D*. Parameter *A* controls how fast the heart can adapt while the terms with HR_{min}, HR_{max}, and *α*, *β* control how difficult it is to reach the maximal heart rate or to rest down to the minimal heart rate.

To learn the parameters of this ODE, previous studies of such models have limited their data collection to controlled laboratory environments with small samples—typically fewer than ten individuals^{14,15,16,17}. Workout intensity data *I* was typically the power output or cadence from a stationary bicycle, and the drive function *f*(*I*) was modeled using low-order polynomials^{14,15}. Here, we extend the ODE method to large-scale uncontrolled environments and use it to model workout data from wearable devices; we do this by learning some parts of the ODE as neural networks and some parameters of the ODE as user-specific variables. Furthermore, we augment the ODE to incorporate environmental factors, such as temperature and fatigue.

### Modeling heart rate in uncontrolled environments

In uncontrolled environments, accurately measuring the exercise intensity *I* becomes challenging. Instead, wearables rely on sensors like GPS and pedometers to measure proxy variables for activity intensity (such as speed, elevation, and number of steps). However, the relationship between *I* and *f*(*I*) becomes unclear when dealing with such generally noisy measurements. To tackle this challenge, we use a flexible neural network to model and learn the drive function *f*.

Uncontrolled environments might also have external factors that can affect heart rate. For instance, higher temperatures induce a higher oxygen demand^{21}. In contrast to controlled environments, workouts in uncontrolled settings may also vary in duration, leading to potential changes in heart rate dynamics over time due to fatigue. Hence, we refine the heart rate equations by modeling how they are affected by the weather *W*, which includes temperature and humidity measurements at the time of workout, and add the effect of fatigue incurred over time *t* during the workout. We parameterize these effects by neural networks *g*(*W*) and *h*(*t*), respectively, and we incorporate them into the demand equation. The term (*f*(*I*(*t*)) − *D*(*t*)) becomes (*f*(*I*(*t*)) ⋅ ** g(W) ⋅ h(t)** −

*D*(

*t*)). For instance,

*g*(

*W*) > 1 indicates an increase in oxygen demand for the weather

*W*.

### Learning a personalized large-scale heart rate model

Each individual possesses their own set of personalized parameters (A, B, *α*, *β*, as well as the drive function *f* and HR_{min}, HR_{max}) that capture their unique heart rate dynamics in response to exercise. Inferring these parameters for each subject and understanding how they might evolve over time can reveal important insight into their health status^{17}. However, learning these parameters for every subject and each new workout is computationally expensive. Instead of directly learning a set of parameters for each subject, we assume that an individual’s health state at a given time can be represented by a low-dimensional latent vector \(z\in {{\mathbb{R}}}^{\ell }\). Then, we turn each ODE parameter into a function of this health representation. For instance, the parameter *α* becomes *α*(*z*), and *f*(*I*) becomes *f*(*z*, *I*). All these “functions of *z*” are parametrized as neural networks, and our goal will be to learn these health representations.

With these changes, the ODE in Equation (1) becomes

$$\left\{\begin{array}{ll}\dot{D}(t)\;\;\;\;=\,B{{{\boldsymbol{(z)}}}}\cdot \left(f({{{\boldsymbol{z,}}}}\,I(t))\cdot {{{\boldsymbol{g(W)}}}}\cdot {{{\boldsymbol{h(t)}}}}-D(t)\right),\\ \dot{{{{\rm{HR}}}}}(t)\;\;=\,A{{{\boldsymbol{(z)}}}}\cdot {\left({{{\rm{HR}}}}(t)-{{{{\rm{HR}}}}}_{{\mathrm{min}}}{{{\boldsymbol{(z)}}}}\right)}^{\alpha {{{\boldsymbol{(z)}}}}}\cdot {\left({{{{\rm{HR}}}}}_{{\mathrm{max}}}{{{\boldsymbol{(z)}}}}-{{{\rm{HR}}}}(t)\right)}^{\beta {{{\boldsymbol{(z)}}}}}\cdot \left(D(t)-{{{\rm{HR}}}}(t)\right),\\ {{{\rm{HR}}}}(0)\;\,=\,{{{{\rm{HR}}}}}_{0}{{{\boldsymbol{(z)}}}}\\ D(0)\;\;\;\,=\,{D}_{0}{{{\boldsymbol{(z)}}}}.\end{array}\right.$$

(2)

Fitting our proposed ODE model to large-scale wearable data involves learning both the global shared neural networks and the subject-time specific health representation *z*. In order to efficiently learn the health representation *z*, we finally introduce one last component inspired by deep learning methodologies. We posit that at any time, the subject’s workout history up to this time contains all of the information to characterize *z*. To model the complex interactions that define *z* based on this history, we utilize a convolutional neural network (CNN) encoder architecture. This allows us to learn the subject-specific health representation as a function of their past workout data. Hence, we exchangeably refer to *z* as the *health representation* or the *history embedding*. Figure 2 summarizes the full pipeline of our method: at any given time and for any subject, the model takes the workout history of this subject up to this time and feeds it to an encoder function to obtain a health representation *z*. Subsequently, the representation *z* is transformed into ODE parameters that are used to solve the ODE for new incoming workouts. Training the model end-to-end is done with standard gradient descent to identify the best neural network weights that best predict workout heart rate sequences.

For training and evaluation, we divided the data into a training set and a testing set. The training set comprises the first 80% of workouts for each subject, while the remaining 20% of workouts form the test set and are held out during training. We selected a few hyperparameters using the best training loss. Additional model details and a description of the implementation, neural network hyperparameters, and training procedure using ODE solvers can be found in Supplementary Methods.