Time dependence of heart rate variability during treadmill running

To investigate the time dependence of the heart rate variability (HRV) during treadmill running, a feedback control loop was implemented to eliminate the potentially confounding influence of cardiovascular drift. Without cardiovascular drift, observed changes in HRV can be directly attributed to time only and not to drift-related increases in heart rate. To quantify the time-dependence of HRV, standard HRV metrics for two consecutive windows of equal duration (12.5 min) were computed and compared. Eight participants were included. The outcome measures showed an overall tendency to decrease over time. Seven of the 10 HRV metrics were significantly lower (p<0.05); three HRV metrics showed moderate evidence of decrease over time, viz. average control power (p = 0.053), very-low frequency power (VLF) of the RR-signal (p = 0.072) and low frequency power (LF) of the RR-signal (p = 0.12). Taken together, these results provide evidence of a decrease in HRV over time during treadmill running; the employment of feedback control of heart rate is important as cardiovascular drift was eliminated. Further work is required to optimize the experimental design and to use a larger sample size to improve the statistical power of the results.


Introduction
Heart rate variability (HRV) is an important physiological indicator related to interactions between the sympathetic and parasympathetic divisions of the autonomic nervous system (Acharya et al., 2005;Kim et al., 2018). Decreases in HRV can be used as a strong predictor of cardiac and arrhythmic mortality (Stein et al., 2005;Sztajzel, 2004). Further, HRV was shown to be able to predict strokes (Lees et al., 2018). With short term (0.5 to 5 minutes) or long term (24 hour) recordings (Bourdillon et al., 2017;Shaffer & Ginsberg, 2017), HRV has become an established, noninvasive method to provide some insight into the activity of the autonomic nervous system (Sztajzel, 2004). Along with the medical field, the field of sports and exercise sciences has shown increased interest in time and frequency domain indices of HRV in recent years (Dong, 2016).
HRV is defined as the variation of the time interval between two successive heartbeats (Shaffer & Ginsberg, 2017). To quantify HRV, the RR-interval is measured: the RR-interval is the time between two detected heart beats, calculated for every R-wave of a QRSevent (Malik et al., 1996). If unreliable RR-intervals (cased by undetected beats, artefacts or arrhythmia) are excluded in post-processing, the signal is commonly referred to as NN-interval (normal-to-normal) (Shaffer & CONTACT Lars Brockmann lars.brockmann@bfh.ch Ginsberg,2017). The exclusion of abnormalities is a crucial step in a sophisticated HRV analysis (Choi & Shin, 2018). Standards for HRV outcome analysis employ time-and frequency-domain measures (Ishaque et al., 2021;Malik et al., 1996). Commonly used time-domain measures are the standard deviation of the NN-intervals (SDNN) and the root mean square of successive RR interval differences (RMSSD). For the analysis of HRV in the frequency domain, there are usually four distinct frequency bands to consider (Malik et al., 1996): • Ultra-low frequency (ULF), with f < 0.003 Hz; • Very-low frequency (VLF), with 0.003 ≤ f < 0.04 Hz; • Low frequency (LF), with 0.04 ≤ f < 0.15 Hz; • High frequency (HF), with 0.15 ≤ f ≤ 0.4 Hz.
The frequency-domain measures estimate the mean power in those four frequency bands (Pagani et al., 1984;Shaffer & Ginsberg, 2017). The power is generally calculated by integrating the spectral density estimate.
The focus of the present work is on changes in heart rate over time during moderate-intensity treadmill running. A recent review systematically analysed evidence relating to changes in HRV with respect to intensity, duration and modality (Michael et al., 2017). The two principal conclusions of the review were: (i) the intensity of exercise is the main factor affecting HRV, where a substantial reduction in HRV occurs as intensity increases; (ii) HRV appears to decrease over time, but only at relatively low intensity and when cardiovascular drift is present (this is the tendency for a progressive decline in stroke volume over time coupled with an increase in heart rate). A limitation of the studies included in the review is that only short-duration exercise sessions were studied, meaning that only LF and HF components of HRV power could be considered.
To overcome these limitations, previous work considered exercise durations sufficient to allow analysis of all classical HRV frequency bands, including VLF and ULF (Hunt & Saengsuwan, 2018): consistent with previous reports (Michael et al., 2017), it was found that HRV decreased with increasing exercise intensity and also that HRV decreased over time. However, in relation to time dependency, it was concluded that 'it remains to clarify whether these changes are due to time itself or to increases in HR related to cardiovascular drift' (Hunt & Saengsuwan, 2018).
A novelty of the approach taken in the present work, specifically introduced to overcome the drift-related limitation identified above, is the use of feedback control to keep heart rate constant (Hunt & Fankhauser, 2016;Hunt & Gerber, 2017), thus eliminating the potentially confounding effect of cardiovascular drift. Having the ability to keep HR constant over the entire duration of an exercise bout enables the unobstructed analysis of HRV's time-dependency, completely independent of HR based influences. The idea of implementing HR control to stabilize the HR for an HRV time-dependence analysis is new and has not been done in any previous research.
The main challenge of this work is to accurately control heart rate in the face of broad-spectrum disturbances arising from the natural, physiological phenomenon of heart rate variability (HRV). But this has to be done in a way that does not excite the control signal (treadmill speed) in a way that is uncomfortable for the runner. The novel inputsensitivity-shaping approach to feedback design that is applied in this work allows the control system design engineer to address these challenges directly. Because classical knowledge on the spectral characteristics of HRV exists, the frequency-domain is the most appropriate framework for feedback design.
Using this novel approach, the work generated new knowledge on the characteristics and dynamics of HRV. This is especially valuable for ultra-low and very-low frequency components, as the understanding is presently still limited. New knowledge represents an important scientific contribution to the field of HR physiology. Furthermore, relevant findings of HRV characteristics can be used to optimize the design of automatic heart rate control systems.
The aim of this work was to investigate the time dependency of HRV during treadmill exercise, under the condition of feedback-controlled heart rate.

Methods
Eight participants were included in the study (Table 1). Of those eight participants, one was female and 7 were male. Inclusion criteria were: male or female, age between 18 and 65 years and physically healthy (exercise 3 × 30 min/wk). Exclusion criteria were: known cardiovascular, pulmonary or musculoskeletal problems that might have interfered with treadmill exercise.

Feedback control
For the control design, a classical feedback structure with a strictly proper linear nominal plant P o and a controller C was employed ( Figure 1). A linear controller design was prioritized over a non-linear one, as to date, there is a lack of evidence that non-linear controllers bring significant performance benefits (Asheghan & Míguez, 2016;Verrelli et al., 2021). Full details of plant model identification and the feedback control design approach can be found elsewhere (Hunt & Fankhauser, 2016), with a brief summary given in the following.
The variable y is heart rate, the control signal u is the speed command sent to the treadmill and the disturbance variable d represents the influence of HRV. For the reference signal r, a constant target heart rate at the border of moderate and vigorous exercise intensity levels, Because of the employment of feedback control, the frequency-domain characteristics of the control loop must be taken into account in relation to the HRV analysis. The principal closed-loop transfer functions to be considered are the sensitivity function S o , which is the transfer function from disturbance d to controlled output y (d → y); the input sensitivity function U o (d, r → u); and the complementary sensitivity function (2). (Here, the reference heart rate r was constant, therefore r and T o are not discussed further in the sequel.) In the experiments, two similar controllers were employed, denoted as C 1 (s) and C 2 (s). These controllers were designed and tested in a previous study that focused on closed-loop feedback control of heart rate (Wang & Hunt, 2021). Because data were available with a constant heart rate target for both of these controllers, both sets of data, with C 1 and C 2 (see Equation (3), with steady-state gains k 1 and k 2 , time constants τ 1 , τ 21 and τ 22 and the feedback design parameter p, which is the bandwidth of U o ), were employed for the analysis in the present work. For these controllers, all sensitivity functions were computed ( Figure 2); it can be seen that the frequency responses with C 1 and C 2 are very similar, and, as a consequence, there was no discernible difference in their performance.
As noted in Equation (2), the sensitivity function and input sensitivity function, respectively, relate the HRV input signal d to the heart rate y and to the control signal The four frequency bands for the heart rate variability analysis are ultra-low frequency (ULF), very low frequency (VLF), low frequency (LF) and high frequency (HF). The red dots mark the respective −3 dB bandwidths. C 1 and C 2 refer to the two controllers employed.
(treadmill speed command) u. The sensitivity function S o governs the influence of the HRV disturbance d on the controlled heart rate y (S o (s): d → y). With the given highpass behaviour of S o (Figure 2), it is apparent that signals with a frequency lower than the −3 dB bandwidth at around 0.0016 Hz, which is close to the border between the ULF and VLF bands, will be attenuated by the compensating action of the control loop. This means that limited ULF activity will be present in the recorded heart rate. A further relevant relationship is described by the input sensitivity function (U o (s): d → u). The low pass characteristic of this function attenuates frequencies higher than the U o bandwidth of p = 0.01 Hz, which is in the middle of the VLF band, thus reducing VLF, LF and HF components of HRV in the signal u.
In view of these frequency-domain characteristics of the control loop, the heart rate signal y can be used primarily to analyse VLF, LF and HF components of HRV (high-pass character of S o ), while the control signal u contains information on the ULF band and, in part, the VLF band (low-pass character of U o ). Furthermore, the information content of u in the LF and HF bands is limited by the data sampling rate (see below, Section 2.4.2).

Experimental design and test procedures
Each participant took part in two tests, one using controller C 1 and one with C 2 , and with each test on a separate day. The order of presentation of C 1 and C 2 , i.e. C 1 then C 2 vs. C 2 then C 1 , was counterbalanced between each participant. Each test had the following procedure: a 15 min warm up at moderate intensity followed by a 10 min break; there was then a formal measurement phase of 30 min, with heart rate controlled to a constant level of HR ref (see Equation (1)); when the 30 min test phase was over a 10 min cool down was initiated. For the HRV analysis, two time windows of duration 12.5 min were defined: w 1 was set to the interval from 300 s to 1050 s and w 2 was from 1050 s to 1800 s; the transient phase in the first five minutes of each test was excluded from the analysis. The different phases of the tests can be seen in a representative data record from one test with one participant (Figure 3).
With each of the eight participants performing two tests, 16 pairs of data samples were obtained for the statistical analysis (i.e. n = 16).

Equipment and data collection
A Simulink Desktop Real-Time application (The Math-Works, Inc., USA) was developed to host the HR feedback control algorithm. The application was deployed on a computer, connected to a treadmill (model Venus, h/p/cosmos Sports & Medical GmbH, Germany) and to a HR receiver module (Heart Rate Monitor Interface, Sparkfun Electronics, USA). Bidirectional communication between the computer and the treadmill was established using the coscom v3 protocol (h/p/cosmos Sports & Medical GmbH, Germany). The HR receiver module was connected to the Matlab/Simulink-based controller using a serial interface with a sampling rate of 1 Hz. The Heart rate was obtained using a chest belt sensor (model T34, Polar Electro Oy, Finland). Since raw RR intervals were not recorded directly, they were reconstructed using the HR signal as follows: The controller was set to run with a sample period of T s = 5 s, thus the HR data was downsampled by averaging five individual values over each sample interval.

Outcome measures
In order to investigate the dependence of HRV on time, time-and frequency-domain outcomes were calculated for windows w 1 and w 2 and the outcomes for these two windows were compared.

Analysis in the time-domain
As the standard deviation of the NN intervals SDNN is closely related to the root-mean-square tracking error RMSE, the RMSE value was used as a proxy for the SDNN HRV metric (Hunt et al., 2019a(Hunt et al., , 2019b, as follows: where HR ref is the constant heart rate reference, Equation (1), and N is the number of samples on the evaluation interval.
In a similar vein, average control signal power P ∇u can be used as a proxy for the HRV metric RMSSD (Hunt et al., 2019a(Hunt et al., , 2019b, viz. ( 6 )

Analysis in the frequency domain
The Lomb-Scargle method for spectral density estimation was used to analyse both the RR data and the control signal (treadmill speed command u). No artefact detection and replacement was needed due to sufficiently high signal quality. For each window (w 1 and w 2 ), the average power in the ULF, VLF, LF, HF frequency bands and in the total frequency range was computed from RR data by integrating the power spectral density estimate obtained from the Lomb-Scargle analysis.
Frequency-domain outcomes were also obtained, and in a similar way, using the control signal u. However, not all frequency bands were included in the analysis of u due to the frequency response of the input sensitivity function U o , as discussed above (Section 2.1), and due to the data sampling rate: because of the applied controller sample period of T s = 5 s, the power density estimation could only be evaluated up to the Nyquist frequency of 0.1 Hz (i.e. half of the sample frequency), which lies within the LF band. Thus, LF and HF components were not included in the analysis of speed signal u.

Statistical analysis
All data except RMSE were log-transformed prior to statistical analysis: since all outcomes, with the exception of RMSE, are derived from squared functions, the data can be expected to follow a log-normal distribution (Hunt & Saengsuwan, 2018). Normality of the sample differences was checked using a Lilliefors test: all data except the average power spectral density estimation in the HF band for the RR time series followed a normal distribution. For the normally-distributed samples a paired t-test was used to analyse the mean differences of the windows w 1 and w 2 and for the non-normal sample (average HF power of the RR time series) a non-parametric paired Wilcoxon signed-rank test was used. The hypothesis of this study was that HRV reduces with time ; therefore, one-tailed tests were employed. The significance level for all tests was set to 5% (α = 0.05).
The statistical analysis was mainly performed using the Statistics and Machine Learning Toolbox in Matlab (The MathWorks, Inc., USA), except that, for the calculation of the confidence intervals for the Wilcoxon test, the statistical computing tool R (R Foundation for Statistical Computing, Austria) was used.

Results
For seven of the ten HRV outcome measures, the values for the second window w 2 were significantly lower than for the first window w 1 (p < 0.05; Table 2). For the remaining three measures, there was moderate evidence of reduction (p = 0.053, 0.072 and 0.012 for P ∇u , VLF power RR and LF power RR, respectively). Taken together, these results point to a reduction in HRV over time.
The accuracy of heart-rate control achieved by the feedback can be quantified using the RMSE, Equation (5), for the total measurement window, i.e. for w 1 and w 2 combined: the average RMSE of the total measurement window of all measurements was 2.04 bpm ± 0.53 bpm (mean ± standard deviation), with a range from 1.25 bpm to 2.99 bpm. The average control signal power P ∇u was 2.33 × 10 −4 (m/s) 2 ±1.10 × 10 −4 (m/s) 2 , range from 0.87 × 10 −4 (m/s) 2 to 5.29 × 10 −4 (m/s) 2 . The average treadmill speed was 2.28 m/s ± 0.61 m/s, ranging from 1.57 m/s to 3.40 m/s. A representative data record from one test with one participant is provided (Figure 3).

Discussion
The aim of this work was to investigate time dependency of HRV during treadmill exercise, under the condition of feedback-controlled heart rate. By comparing two consecutive windows of equal duration, significant decreases in the HRV metrics RMSE, ULF power RR, HF power RR, total power RR, ULF power speed, VLF power speed and total power speed were found (p < 0.05). For P ∇u , VLF power RR and LF power speed, the statistical analysis found only moderate evidence of decreases (p = 0.053, p = 0.072 and p = 0.12, respectively). Taken together, these results provide evidence of a decrease in HRV over time during treadmill running.
For the first time, feedback control of heart rate was employed in the analysis of time dependency of HRV. The employment of feedback control of heart rate to eliminate cardiovascular drift has been proven technically feasible. To keep heart rate constant, treadmill speed Table 2. Outcome measures for windows w 1 and w 2 ; n = 16; w 1 : first window; w 2 : second window; SD: standard deviation; MD: mean difference of w 1 -w 2 ; 95 % CI: confidence interval for the mean difference; p-values: paired single-sided t-tests, except for HF power of the RR signal -Wilcoxon signed-rank test; RMSE: root-mean-square error; all outcomes with the exception of RMSE are log transformed and are dimensionless. was automatically and continuously adapted by the feedback and by this means cardiovascular drift was eliminated. This contribution is important because it is known that HRV decreases with increasing heart rate (Michael et al., 2017), which can confound the investigation of time effects only. Thus, the observed decreases in HRV can be directly explained by the influence of time itself. This study has several limitations that can be addressed in future work. Firstly, some aspects of the feedback control design can be modified to improve information content in the measured signals. The lowpass nature of the input sensitivity function U o meant that the speed signal u contained components primarily in the ULF and VLF bands, while the LF and HF content was strongly attenuated. Furthermore, as described above (Section 2.4.2), the Nyquist frequency of 0.1 Hz lies within the LF range and gave a hard upper limit for the analysis of u. By employing alternative design parameters as detailed elsewhere (Hunt & Liu, 2018), it is possible to make the frequency response of U o constant over all frequencies; if, in addition, a higher sampling rate is used, it would then be possible to analyse all HRV bands using the speed signal u.
Secondly, to improve the resolution of the RR signal, raw RR intervals should be recorded directly rather than by reconstruction from the heart rate signal. This modification would likely improve analysis of LF and HF components of the RR signal.
Thirdly, the statistical power of the results can be strengthened in future studies by using a larger sample size; the effect sizes and dispersions observed in the present work can be used in a statistical power analysis to estimate the appropriate sample sizes.
Reflecting on the exercise intensity of this study, it is difficult to evaluate whether the ventilatory frequency stayed largely the same or increased throughout the duration of the exercise. An increasing ventilatory frequency might directly impact the HRV. With a moderate training intensity, however, ventilation would be expected to remain stable over time. Measuring the ventilation rate would resolve any uncertainty by aiming to ensure that the ventilatory frequency remains constant. Another possibility to reduce the influence of ventilation on HRV is to lower the exercise intensity to increase the likelihood of a stable ventilatory frequency.
The HF band used in this work has an upper limit of 0.4 Hz. To include the possibility of higher-frequency components of HRV, for example, related to heart rate and ventilation at heavy exercise intensities in further studies, an additional band could be considered above 0.4 Hz.

Conclusions
The results of this study provide evidence of a decrease in HRV over time during treadmill running; the employment of feedback control of heart rate is important as cardiovascular drift was eliminated. Further work is required to optimize the experimental design and to use a larger sample size to improve the statistical power of the results.

Ethics approval and consent to participate
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. 2019-02184). All participants provided written, informed consent.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This study was funded by the Swiss National Science Foundation as part of the project 'Heart Rate Variability, Dynamics and Control During Exercise' (Ref. 320030-185351).

Data availability statement
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.