Lars Brockmann,Kenneth J. Hunt
Heart rate variability changes with respect to time and exercise intensity during heart-rate-controlled steady-state treadmill running
Biomedical engineering,Biomarkers,Cardiology,Control theory
The aim of this work was to investigate the time and exercise intensity dependence of heart rate variability (HRV). Time-dependent, cardiovascular-drift-related increases in heart rate (HR) were inhibited by enforcing a constant heart rate throughout the exercise with a feedback control system. Thirty-two healthy adults performed HR-stabilised treadmill running exercise at two distinct exercise intensity levels. Standard time and frequency domain HRV metrics were computed and served as outcomes. Significant decreases were detected in 8 of the 14 outcomes for the time dependence analysis and in 6 of the 7 outcomes for the exercise intensity dependence analysis (excluding the experimental speed-signal frequency analysis). Furthermore, metrics that have been reported to reach an intensity-dependent near-zero minimum rapidly (usually at moderate intensity) were found to be near constant over time and only barely decreased with intensity. Taken together, these results highlight that HRV generally decreases with time and with exercise intensity. The intensity-related reductions were found to be greater in value and significance compared to the time-related reductions. Additionally, the results indicate that decreases in HRV metrics with time or exercise intensity are only detectable as long as their metric-specific near-zero minimum has not yet been reached.
Heart rate variability changes with respect to time and exercise intensity during heart-rate-controlled steady-state treadmill running

The aim of this work was to investigate the time and exercise intensity dependence of heart rate variability (HRV). Time-dependent, cardiovascular-drift-related increases in heart rate (HR) were inhibited by enforcing a constant heart rate throughout the exercise with a feedback control system. Thirty-two healthy adults performed HR-stabilised treadmill running exercise at two distinct exercise intensity levels. Standard time and frequency domain HRV metrics were computed and served as outcomes. Significant decreases were detected in 8 of the 14 outcomes for the time dependence analysis and in 6 of the 7 outcomes for the exercise intensity dependence analysis (excluding the experimental speed-signal frequency analysis). Furthermore, metrics that have been reported to reach an intensity-dependent near-zero minimum rapidly (usually at moderate intensity) were found to be near constant over time and only barely decreased with intensity. Taken together, these results highlight that HRV generally decreases with time and with exercise intensity. The intensity-related reductions were found to be greater in value and significance compared to the time-related reductions. Additionally, the results indicate that decreases in HRV metrics with time or exercise intensity are only detectable as long as their metric-specific near-zero minimum has not yet been reached.

Heart rate variability (HRV) is the variation in the time interval between consecutive heartbeats (RR intervals—the time elapsed between two successive R waves of the electrocardiogram QRS complex). It is an important physiological indicator related to interactions between the sympathetic and parasympathetic divisions of the autonomic nervous system (ANS). HRV has become an established, non-invasive indicator used across various fields. A range of methods for measurement, signal analysis, and interpretation of HRV have been developed and established as formal standards, comprising time domain, frequency domain, and joint time-frequency domain approaches. To accurately assess the function of the ANS, time and frequency domain HRV analysis methods are most commonly used. In the frequency domain, four distinct frequency bands have classically been investigated: ultra-low frequency (ULF), with Hz; very-low frequency (VLF), with Hz; low frequency (LF), with Hz; and high frequency (HF), with Hz. A recent review systematically analysed, among other things, cardiac autonomic responses during exercise using HRV, with a focus on intensity, duration, and modality. It was concluded that the intensity of exercise is the main factor affecting HRV, where a substantial reduction in HRV occurs as intensity increases, and that HRV appears to decrease over time for relatively low-intensity exercise. A distinction could not be made as to whether the decrease in HRV over time was caused by the exercise duration (i.e. time) or by the cardiovascular-drift-related elevated heart rate (HR). A limitation of the studies included in the review is that only short-duration exercise sessions were studied, implying that only low-frequency (LF) and high-frequency (HF) components of HRV power could be considered. To overcome these limitations, our previous work considered exercise durations that allow analysis of all classical HRV frequency bands including very-low (VLF) and ultra-low frequency (ULF). Consistent with previous reports, HRV was found to decrease with increasing exercise intensity, and HRV was observed to decrease over time. But the confounding effects of cardiovascular drift on the time-dependent HRV decrease were still at play, and the study concluded that it remains to be clarified whether these changes are due to time or due to increases in HR related to cardiovascular drift. In the work proposed here, to overcome the drift-related limitations, a feedback control system was employed to keep the heart rate constant throughout the duration of treadmill exercise. The technical feasibility of using HR control for the investigation of time-related changes in HRV was proven previously in a pilot study with eight participants; while the outcome measures showed an overall tendency to decrease over time, the decrease was significant () for only seven of the 10 HRV metrics, pointing to the study being statistically underpowered. Further limitations identified in the pilot study, and which are eliminated in the present work, included the employment of a feedback design with low-pass characteristics, thus attenuating information content in the LF and HF bands, and the reconstruction of RR intervals from the HR signal, rather than from direct recordings of raw RR intervals. By forcing the heart rate to stay constant for the duration of the treadmill exercise, confounding HR-related influences can be removed, and the unobstructed HRV analysis becomes viable. A key challenge with this new HR control approach for HRV analysis is the effect the control loop’s compensating actions can have on specific frequency ranges, potentially affecting frequency-related HRV metrics. Following a suggestion made in our pilot work, a distinct feedback control structure was implemented here to address this limitation by using an input-sensitivity-shaping approach to obtain a uniformly flat frequency response across all frequencies of the treadmill speed signal. Therefore, the treadmill speed signal might be able to act, instead of or as well as the recorded RR intervals, as an indirect source for the frequency-domain HRV analysis, where frequency-specific control loop compensation effects are eliminated. The viability of using speed in this way is explored in the sequel. The fundamental idea of implementing feedback control of HR to stabilise HR for time-dependency analysis of HRV, as well as the concept of using the treadmill speed signal as a proxy for the RR interval signal to perform an indirect frequency-domain HRV analysis, is new. Aside from our pilot study, this has not been investigated in any previous research. The aim of this work was to investigate the time and exercise intensity dependence of HRV during steady-state treadmill running while using feedback control to prevent HR drift. Based on the previous findings of Michael et al. and Hunt et al., we hypothesised that HRV can be expected to decrease with increasing exercise intensity and throughout the duration of the exercise.

Thirty-two healthy, regularly exercising (three times per week and at least 30 min per session) adults were included in the study (Table 1; details of an a priori sample size estimate are given later, “Sample size estimate”). Recruitment was carried out by convenience sampling within the Department of Engineering and Information Technology of Bern University of Applied Sciences in Burgdorf and the University of Bern. Of the 32 participants, 29 were male and 3 were female. Smokers and persons with prior history of cardiovascular or respiratory disease or current musculoskeletal complaints or injuries were excluded. Before each test, participants were required to avoid strenuous exercise (24 h), caffeine (12 h), and heavy meals (4 h). The study was approved by the local ethics committee (Ethics Committee of the Swiss Canton of Bern, Ref. 2021-00889), and the participants provided written informed consent prior to participation.

The control structure was set up as a generic negative feedback control system (Fig. 1) with feedback compensator C(s) and nominal plant .

The nominal plant representing the dynamic response between the treadmill speed u and the heart rate y was defined as a strictly proper transfer function , constrained in the following to be of first order,where k (units beats-per-minute/(m/s)) is the steady-state gain and (units s) is the time constant; and are the degrees of polynomials and . It is known that the plant parameters k and vary widely between different people, yet it has been clearly demonstrated that a constant-coefficient, linear, time-invariant compensator based on an average nominal model can deliver stable and accurate HR control performance when applied to a wide range of participants, that is to say, robust control is obtained as a consequence of the fundamental ability of feedback to reduce plant uncertainty. The nominal plant parameters employed here were set to average values obtained from a total of 73 participants from two separate model identification studies, namely and . These values were derived by combining parameters obtained empirically in two previous studies, more specifically with n = 25 and with n = 48; the respective sample sizes were taken into account. The linear feedback compensator, generally described in rational form as , was intentionally constructed as a merely proper first-order transfer function with an integral term, viz.where and are the degrees of polynomials G and H. As detailed below, this combination of nominal plant Eq. (1) and compensator Eq. (2) results in a closed-loop characteristic polynomial of degree two and, by virtue of the compensator’s two free parameter and , allows arbitrary placement of the two closed-loop poles in the complex plane. The choice of a proper transfer-function for C is to ensure that the input sensitivity function (Eq. (5), below) is also proper and thus remains finite over all frequencies, in line with our design goal to make constant for all frequencies (a strictly proper C would make strictly proper, whence would roll off towards zero as frequency increases). The closed-loop sensitivity function , complementary sensitivity function and the input sensitivity function can be deduced from the control structure. describes the transfer function from disturbance d to the controlled output y (), from the reference signal r to the controlled output y (), and from the reference signal r and the disturbance d to the treadmill speed signal u ( and ): The closed-loop characteristic polynomial can be identified as The feedback design goal is to shape the input-sensitivity magnitude to be constant across all frequencies. The rationale for this goal is that the treadmill speed command signal u, which is linked to HRV disturbance term d through , might potentially be used for frequency-domain HRV analysis across the whole frequency range. From Eq. (5), , and we proceed by taking a cancellation approach to simplify this down to a constant value. As a first step, all plant poles are cancelled by constraining the compensator numerator polynomial G to include the known factor and a remaining, unknown factor by writing , which results in Since the plant is assumed strictly proper, , and the compensator proper, , it follows that the order of the characteristic polynomial, Eq. (6), is . Since a unique solution of Eq. (6) for G and H requires that be equal to the number of free parameters in G and H, namely , it follows that (giving also , since ). This in turn implies that, since by construction , the polynomial must be of degree zero, i.e. is a constant, which we denote . This allows to be written, using Eq. (7), as The next step is to place the remaining closed-loop poles (the roots of ) at the open-loop pole locations by setting Substituting in Eq. (8), this results finally, and as desired, in a constant ,whereby compensator synthesis amounts to solving Eq. (9) for and H. Thus far, the derivation was based on the general plant and compensator , but we now specialise the solution to the first order case defined in Eqs. (1) and (2), i.e. using , , and . With these values, the design equation (9) becomesgiving the solution . It follows from Eq. (10) that The various cancellations involved in the derivation also lead to greatly simplified forms for the sensitivity and complementary sensitivity functions in Eqs. (3) and (4), i.e.and These expressions show that both and have the same bandwidth as the open-loop plant (cf. Eq. (1)), and, indeed, that the closed-loop transfer function is the same as the open loop (up to the scaling factor 1/k): these observations point to the “neutrality” of a compensation strategy that achieves a constant input sensitivity function of magnitude equal to the inverse of the steady-state plant gain. The compensator parameters in can be identified by noting to give and , whereforeand it is seen that the compensator parameters depend only on the plant parameters. Using the nominal values and , the specific compensator used in this study is calculated as The sensitivity functions , and for this compensator with the nominal plant parameters can be computed using Eqs. (12)–(14) and are plotted in Fig. 2. In particular, it can be seen that (blue line) maintains a constant value of 1/k over all frequencies ( and dB) and that the bandwidth of both and is the same as the open-loop plant bandwidth ( and Hz). We remark, for the purposes of later discussion, that this bandwidth is just below the upper bound of the ULF frequency band which lies at Hz.

Each participant performed two treadmill running tests at different exercise intensities (lower intensity level 1 [] and higher intensity level 2 []), and each with a duration of 35 min. The intensities for the two tests were set for each participant individually. Each test was carried out on a separate day with at least 48 h between tests. The recording of the two tests ( and ) was counterbalanced by randomizing the order of employed intensity in order to avoid any order-of-presentation effects. A 10-min warm up was used to assess the HR level and speed at a given baseline exertion level. Participants were asked to choose a running speed equivalent to an exertion level of 13 (i.e. somewhat hard) on the Borg rating of perceived exertion (RPE) scale. The HR recorded during the warm up () was then averaged and used to calculate the two exercise intensities and . The reference HR for each EIL was calculated as a symmetric deviation from by of the maximal age-related heart rate () as corresponds to the lower exercise intensity level and to the higher one. Both and measurements were performed as feedback control tests, where the heart rate was kept constant throughout the running exercise.

Time dependence analysis was performed by calculating and comparing HRV metrics for two consecutive evaluation windows ( and ) with equal duration of 15 min: window was defined from to and window from to . For the exercise intensity dependence analysis, HRV metrics were computed for the total duration (denoted ) of and then compared to the HRV metrics computed for the total duration () of . RR interval outliers or ectopic heartbeats were identified and removed using an impulse rejection filter with a threshold value of 7 and a window length calculated to span about 25 s.

In the time domain, the root-mean-square of successive differences (RMSSD) and the standard deviation of the normal-to-normal (NN) intervals (SDNN) were evaluated and served as primary endpoints. SDNN and RMSSD are defined as follows:where is the ith recorded interval and is the mean value of all intervals. N is the total number of intervals.

In the frequency domain, the average power contained in the frequency bands ULF, VLF, LF, HF, and the total frequency range (TP, Total Power), was computed for the RR signal and the treadmill speed signal using the Lomb-Scargle method for spectral density estimation.

An a priori statistical power calculation using observed effect size and dispersion from a pilot study with a sample of participants was performed to obtain a sample-size estimate. Calculations used a significance level of and a required statistical power of . The statistical power calculation was based on the total power of the HRV spectrum (TP, total power spectral density estimate). The pilot study found a difference in paired samples for the log-transformed data of (mean ± standard deviation). The power calculation revealed a required sample size of . To cater for the possibility of participants dropping out, a 10% contingency buffer (+3 participants) was added, resulting in a final sample size of .

Based on our hypothesis that HRV can be expected to decrease with increasing exercise intensity and throughout the duration of the exercise, one-sided t-tests with significance level were performed for all HRV metrics (time and frequency domain) to assess time and exercise intensity dependence. For the time dependence analysis, HRV metrics computed for were compared to the metrics for (for both intensity levels and ). For the exercise intensity dependence analysis, HRV metrics for the total duration of were compared to those for the total duration of . All metrics were log-transformed to reduce skew and make the data conform more closely to the normal distribution. A Shapiro–Wilk test was performed to assess normality, and no data was found to be significantly different from normal.

All tests were performed on a computer-controlled treadmill (model Pulsar, h/p/cosmos Sports & Medical GmbH, Germany). Treadmill speed was set by a computer running the HR feedback controller within the Simulink Desktop Real-Time environment (The MathWorks, Inc., USA) and communicating with the treadmill over the coscom v3 interface protocol. Heart rate and raw RR intervals were recorded with a chest-strap-mounted sensor (H10, Polar Electro Oy, Finland) transmitted to an ESP32 development board (Espressif Systems, China) over Bluetooth low energy and sent via serial communication to the control application running on the PC. A Polar V800 wristwatch was employed as a backup method for saving the heart rate and RR intervals. The control application was set up to work with heart rate values transmitted in units of beats per minute (bpm). The V800 and the ESP32 development board saved the more accurate RR intervals, recorded with millisecond resolution, for later analysis.

This research was performed in accordance with the Declaration of Helsinki. The study was reviewed and approved by the Ethics Committee of the Swiss Canton of Bern (Ref. 2021-00889). All participants provided written informed consent.

Sixty-two of the 64 measurements were recorded successfully (one participant dropped out of the study due to muscle pain). Fifty-five of the 62 recorded data sets were able to be used for analysis: three measurements had to be excluded due to insufficient HR control performance, another three due to HR abnormalities and one due to faulty sensor data. For illustration, a sample data set for a single participant is provided (Fig. 3), while a complete set of categorical scatter plots showing all hypothesis testing results is provided for time domain (Fig. 4) and frequency domain (Fig. 5) outcomes. Table 2 lists the corresponding numerical values.

SDNN was found to significantly decrease over time at () and (), and was observed to be lower at a higher intensity (). RMSSD, on the other hand, revealed no significant differences, neither for the intensity comparison () nor the time dependence analysis at (), nor the time dependence analysis at ().

The frequency domain outcomes are separated into results derived from the RR interval analysis and results derived from the speed signal analysis (Fig. 5). For the RR interval frequency domain analysis, all frequency-related HRV metrics (ULF, VLF, LF, HF, and TP) were found to decrease significantly with increasing exercise intensity (; ; ; ; ). The outcomes for the time dependence analysis at showed significant decreases for ULF, VLF, and TP (; ; ), moderate evidence for LF () and little to no evidence for HF (). For the time dependence analysis at , LF, VLF, and TP were observed to decrease significantly (; ; ). ULF showed moderate evidence of decrease (). For HF, no evidence of a decrease was found (). For the speed signal frequency domain analysis, ULF, VLF, and TP were found to decrease over time at both and (; ; for and ; ; for ). A significant difference was found in LF at (), and moderate evidence for a decrease was observed in LF at (). No significant differences in HF at () and moderate evidence at () were identified. The exercise intensity dependence analysis revealed significant differences in LF and no evidence in ULF, VLF, HF, and TP (; ; ; ; ).

This work aimed to investigate the time and exercise intensity dependence of HRV during steady-state treadmill running while using feedback control to prevent HR drift: we hypothesised that HRV can be expected to decrease with increasing exercise intensity and throughout the duration of the exercise. Analysing the time dependence of time and frequency domain HRV metrics, and revealed similar findings. SDNN, as well as ULF, VLF, LF, and TP (all except HF), computed using RR intervals and the treadmill speed signal, in the main showed strong evidence of decreasing in value over time (with moderate evidence for at ; at and at ). RMSSD and HF (RR and speed) did not significantly decrease over time. Combining these results with the tendency that most HRV measures reach an intensity-dependent near-zero minimum, and that HRV metrics believed to be associated with cardiac parasympathetic activity (i.e. RMSSD and HF) have been reported to reach that minimum more rapidly (usually at moderate intensity), gives good reason to infer that HRV decreases with time as long as the intensity-dependent near-zero minimum has not yet been reached. This conclusion is also supported by the findings from our pilot study, where exercise performed at light intensity revealed HF and an RMSSD-related metric to decrease over time ( for RMSSD proxy; for HF). It would seem that with the low exercise intensity, RMSSD and HF had not yet reached a minimum and thus were able to decrease over the observation interval. The intensity-level comparison identified a significant decrease in SDNN and in all frequency-domain HRV metrics (computed for the RR intervals) with increasing exercise intensity. RMSSD did not change, indicating the presence of an intensity-dependent near-zero minimum below . Despite the consistent decrease in frequency-domain HRV metrics computed for the RR intervals, the corresponding outcomes for the speed signal’s frequency-domain HRV metrics differed greatly. VLF and HF stayed approximately constant, while ULF and TP increased. This finding is contrary to expectations considering the connection of the RR interval signal via the controller to the speed signal. With the control loop’s flat input sensitivity function and the consistent control structure across both exercise intensities, intensity-dependent changes in outcomes were expected to be reflected in both the RR intervals and the speed signal. Based on the large increase in ULF power, we suspect this outcome to be affected by an inconsistent plant model. A deviation from the nominal plant can directly affect the actual characteristics of the sensitivity functions. The conspicuous increase in ULF leads us to believe that a decrease in the actual steady-state gain parameter k with increasing exercise intensity might have been the cause as the limit of as frequency tends to 0 is 1/k (see “Feedback control design”). A lower than nominal k would reduce the control loop’s dampening effect on the impact the RR signal can have on the treadmill speed, thus leading to an overall increase in power. However, this argument may not be in agreement with a previous study, where model plant parameters k and were reported not to be significantly dependent on exercise intensity. A second contributory factor might have been the larger speed reduction observed at higher intensities. Throughout a running exercise bout, the treadmill speed typically drifts downwards to compensate for cardiovascular-drift-related HR changes: this downward trend in treadmill speed was observed to be greater at higher intensities. As a consequence, we suggest that the method of using the treadmill speed signal as a proxy for the RR intervals during heart rate stabilised running exercise, in order to perform frequency-domain HRV analysis where frequency-specific control loop compensation effects are eliminated using input-sensitivity-shaping, needs further investigation to be applicable for an intensity dependence analysis. On the other hand, for the time dependence analysis, this method produced results that closely matched the trends found in frequencies primarily unaffected by the control loop (VLF, LF, and HF) of the original RR interval analysis. This suggests that the speed signal could be used as an alternative to detect trends in frequencies primarily affected by the control (namely ULF) without the influence of the control loop’s compensation effects. We acknowledge that the duration of 15 min for the analysis windows and is relatively short for the ULF power estimation. Longer analysis windows are generally more desirable but would place additional demands on the participants performing the running exercise, therefore a balance has to be found.

In summary, feedback control of heart rate was successfully employed to answer the question of whether time-dependent HRV changes occur due to time itself or due to cardiovascular-drift-related heart rate increases. Most HRV metrics were found to decrease with time and with exercise intensity. The exercise-intensity-related reductions were generally found to be greater in value and significance compared to the time-related reductions. HRV metrics that have been reported to reach an intensity-dependent near-zero minimum rapidly (usually at moderate intensity) were found to be near constant over time and only barely decreased with intensity, indicating that decreases in HRV metrics with time or exercise intensity are only detectable as long as their metric-specific near-zero minimum has not yet been reached.