 Research
 Open Access
 Published:
Identification of heart rate dynamics during treadmill exercise: comparison of first and secondorder models
BioMedical Engineering OnLine volume 20, Article number: 37 (2021)
Abstract
Background
Characterisation of heart rate (HR) dynamics and their dependence on exercise intensity provides a basis for feedback design of automatic HR control systems. This work aimed to investigate whether the secondorder models with separate Phase I and Phase II components of HR response can achieve better fitting performance compared to the firstorder models that do not delineate the two phases.
Methods
Eleven participants each performed two openloop identification tests while running at moderatetovigorous intensity on a treadmill. Treadmill speed was changed as a pseudorandom binary sequence (PRBS) to excite both the Phase I and Phase II components. A counterbalanced crossvalidation approach was implemented for model parameter estimation and validation.
Results
Comparison of validation outcomes for 22 pairs of first and secondorder models showed that rootmeansquare error (RMSE) was significantly lower and fit (normalised RMSE) significantly higher for the secondorder models: RMSE was 2.07 bpm ± 0.36 bpm vs. 2.27 bpm ± 0.36 bpm (bpm = beats per min), second order vs. first order, with \(p = 2.8 \times 10^{10}\); fit was \(54.5\% \pm 5.2\)% vs. \(50.2\% \pm 4.8\)%, \(p = 6.8 \times 10^{10}\).
Conclusion
Secondorder models give significantly better goodnessoffit than firstorder models, likely due to the inclusion of both Phase I and Phase II components of heart rate response. Future work should investigate alternative parameterisations of the PRBS excitation, and whether feedback controllers calculated using secondorder models give better performance than those based on firstorder models.
Background
Characterisation of heart rate (HR) dynamics with respect to changes in exercise intensity provides models that can be used to synthesise control algorithms to maintain target HR levels [1]. The control of HR is important in the design of training protocols that aim both to maintain and to improve cardiorespiratory fitness; this applies to healthy individuals [2] and also in different patient populations [3, 4]. Target heart rate profiles come in various forms such as highintensity interval training (HIIT) that repeats highintensity exercise connected by lowintensity recovery intervals; HIIT has potential to enhance cardiovascular health and fitness when compared to training at constant work rates (systematic reviews: [5, 6]).
Recent work investigated the effect of exercise intensity and time on HR dynamics using firstorder models [7], but it may be beneficial to include higher order effects: based on physiological study, the dynamics of both oxygen uptake and HR responses to changes in exercise intensity are known to have three distinct phases [8]. These are: (i) a Phase I component lasting \(\sim \) 15 s with a relatively smallmagnitude ventilatory response, but where HR can increase by about 50% of its total response [9]; (ii) a Phase II component between around 15 s and 3 min contributing the further increase of cardiopulmonary response; and then, (iii) if the applied exercise intensity exceeds the anaerobic threshold, a Phase III component is prolonged and rises slowly. The three components can each be modelled as single exponentials (firstorder systems) each with their own time delay, gain, and time constant [10]. In addition to these primary dynamic responses, the phenomenon of heart rate variability (HRV) can be added to the model to represent the regulatory activities of the autonomic nervous system; in the context of feedback control of HR, HRV represents a broadspectrum disturbance term [1].
Because it can be challenging to estimate the separate Phase I and II components using data which is noisy, those two phases have often been identified as a combined single exponential model with a time constant termed the mean response time (MRT, [8]), which is effectively the firstorder approach taken in the previous studies that focused on system identification [7] and feedback control [1]. In feedback control, the slow Phase III component can readily be neglected as it is compensated by inclusion of an integrator in the controller. The focus of the present work is therefore the investigation of whether the separate identification of Phase I and II components, i.e., the employment of a secondorder model, can give better model fidelity.
Other recent approaches to HR dynamics identification focused mainly on the modelling of the Phase II and III components of the HR response. Several studies employed a nonlinear statespace model structure comprising two different states (\(x_1\) and \(x_2\)) to separately describe the Phase II and Phase III dynamics [11,12,13,14,15]. Other work used linear timevarying systems to model the slow Phase III dynamic [16, 17]. While inclusion of Phase III may improve overall model fidelity, it will, as noted above, have negligible impact on feedbackcontrol performance as it will be eliminated by the integral action. In contrast, it can be anticipated that separate modelling of the Phase I and II components might lead to better control performance when the model is used as the basis of an analytical, modelbased feedback design.
To this end, this work aimed to investigate whether secondorder models with separate Phase I and Phase II components of HR response can achieve better fitting performance compared to firstorder models that do not delineate the two phases. Here, an input signal of PRBS (pseudorandom binary sequence) form was designed to excite both the Phase I and Phase II components.
Results
To illustrate the procedures of data preprocessing and model validation, an exemplary result from participant P04 is shown (Fig. 1); the raw data for the same participant are shown above below in the section ‘Method’. For this example, the secondorder model \(P_2\) gave better performance than the firstorder model \(P_1\): fit was 51.9% vs. 50.9% (\(P_2\) vs. \(P_1\)) and RMSE was 2.01 bpm vs. 2.05 bpm.
The overall statistical comparison of validation outcomes for the 22 pairs of first and secondorder models showed that RMSE was significantly lower and fit significantly higher for the secondorder models: RMSE was 2.07 bpm ± 0.36 bpm vs. 2.27 bpm ± 0.36 bpm, \(P_2\) vs. \(P_1\), with \(p = 2.8 \times 10^{10}\) (Table 1; Fig. 2a); fit was \(54.5\% \pm 5.2\)% vs. \(50.2\% \pm 4.8\)%, \(p = 6.8 \times 10^{10}\) (Table 1; Fig. 2b). The graphical illustration of overall outcomes (Fig. 2) shows how widely individual samples and their differences are dispersed, together with means and their 95% confidence intervals (CIs). These plots allow visual determination of significant differences, if they exist: whenever there is a significant difference, the value 0 will not be contained within the corresponding CI.
The sample size was estimated a priori by a statistical power calculation that used estimates of expected effect sizes and sample standard deviations, with significance level set to 5% (\(\alpha \) = 0.05) and with a statistical power of 80% (\(1\beta \) = 0.8).
The observed outcomes show large effect sizes (approximately 9% for both outcomes) and extremely low p values (on the order of \(10^{10}\)), thus pointing to a wellpowered statistical analysis. In fact, post hoc statistical power analysis based on observed effect sizes and sample dispersions gives an observed power of 100% for both outcomes.
A graphical illustration of the dispersion of estimated model parameters for the 22 first and 22 secondorder models is provided (Fig. 3). The overall first and secondorder models were obtained by averaging the individual gains and time constants. For the firstorder models, the overall gain was \(k_1\) = 28.57 bpm/(m/s) ± 5.27 bpm/(m/s) (mean ± standard deviation) while the time constant was \(\tau _1\) = 70.56 s ± 16.84 s. For the secondorder models, the overall gain was \(k_2\) = 24.70 bpm/(m/s) ± 5.07 bpm/(m/s) and the overall time constants were \(\tau _{21}\) = 18.60 s ± 7.88 s and \(\tau _{22}\) = 37.95 s ± 16.01 s. This gives the average transfer functions for first and secondorder models as follows:
Discussion
This study aimed to investigate whether secondorder models with separate Phase I and Phase II components of heart rate response can achieve better fitting performance compared to firstorder models that do not delineate the two components. The results clearly demonstrate that secondorder models give significantly better goodnessoffit, in terms of both RMSE and fit (NRMSE): RMSE was on average 0.19 bpm lower and fit 4.3% higher for the secondorder model structure (p values were on the order of \(10^{10}\) in both cases); that these significance levels were achieved with a small sample size of only 11 participants underline the difference.
The approach taken here focused on controlorientated model structures, in the sense that the estimated models would be intended to be used for analytical (modelbased) design of heart rate control systems. For this reason, slow Phase III components in the data were eliminated by detrending prior to parameter estimation. This is consistent with feedbackcontrol scenarios where slowly drifting Phase III variations in heart rate are automatically compensated using integral action in the controller.
A further difference between the methodology employed here and heart rate modelling approaches taken in the physiological literature, [8], is that a nominal operating point was assumed, and small deviations around this point were considered (in this case, the operating point was set at the transition between exercise levels considered to be moderate and vigorous). This is consistent with linear feedback design approaches, which are implicitly based on models that are smallsignal linearisations around an operating point; the purpose of feedback control is indeed to maintain the controlled variable, viz., heart rate, close to a target level.
For these reasons, it is not possible to compare the overall estimated model parameters (gains and time constants, Eqs. (1) and (2) with values given in the physiological literature (e.g., [9, 10]), because, there, responses are usually recorded using large steps from a resting or lowintensity baseline.
A consequence of the controlorientated methodology followed here is that the design of the PRBS input signal becomes important. For nonlinear systems, it is known that the parameters of linear approximations are input dependent [18], which motivates further work to explore the effect of PRBS amplitude and frequency content on model fidelity; in particular, it is important to focus the information content on frequencies around the intended crossover band of the closedloop system [19].
Future work should investigate whether the observed improvement in model fidelity translates into better feedbackcontrol performance, i.e., whether controllers designed on the basis of secondorder models perform better, in some sense, than those designed using firstorder models. Because of the fundamental property of feedback that plant uncertainty (including modelling error) is reduced, the answer to this question will likely not be as clear cut as in the openloop identification case.
Conclusions
Secondorder models give significantly better goodnessoffit than firstorder models, likely due to the inclusion of both Phase I and Phase II components of heart rate response. Future work should investigate alternative parameterisations of the PRBS excitation, and whether feedback controllers calculated using secondorder models give better performance than those based on firstorder models.
Methods
Participants
Eleven healthy participants were recruited (8 males, 3 females) with age 32.5 years ± 12.3 years (mean ± standard deviation), body mass 75.5 kg ± 14.4 kg, and height 179 cm ± 12 cm. For inclusion, each participant was required to be a regular exerciser (30min bouts, 3 times per week) and nonsmoker, and to be free of injury and illness.
Test protocols
To generate separate estimation and validation data sets, each participant took part in two identification tests; there was an interval of at least 48 h between the two tests. Before each test, participants were asked to meet the following requirements: refrain from strenuous activity for 24 h, caffeine for 12 h, avoid large meals for 3 h. Each test session had four phases: a 15 min warm up, a 10 min rest, a 36 min formal measurement, and a 10 min cool down (Fig. 4a).
In the warm up, a feedbackcontrol system was employed to automatically regulate the speed of the treadmill to maintain a constant target HR. The target HR, denoted \(\mathrm{HR}_{\mathrm{ref}}\), was computed individually for each participant and corresponded to the HR at the transition between intensity levels considered to be moderate or vigorous [3], as follows: \(\mathrm{HR}_{\mathrm{ref}} = 0.765\times (220  \mathrm{age})\) [beats/min, bpm] (except for participant P03, for whom the factor 0.7 was used, because 0.765 led to HR remaining in the vigorousintensity regime). The mean speed of the treadmill during the final 2 min of the warm up phase was subsequently used as the midlevel speed, denoted \(v_\text{m}\), for the next phase.
In the formal measurement phase, the speed of the treadmill, denoted v, was designed as a fifthorder PRBS with mean speed \(v_{\rm m}\) and amplitude 0.25 m/s, i.e., \(v = v_{\mathrm m} \pm 0.25~\) m/s (to illustrate, a single original data record is provided; Fig. 4b). Model parameter estimation and validation was performed over a full cycle of the PRBS using an evaluation period from 290 s to 2085 s (Fig. 4); the first 5 min were excluded to eliminate the initial transient. During the cool down phase, the speed of the treadmill was kept constant at \(v = v_{\mathrm{m}}  0.5~\) m/s.
Equipment
All tests were carried out using a treadmill (model Venus, h/p/cosmos Sports & Medical GmbH, Germany) controlled by a PC running realtime Matlab/Simulink (The MathWorks, Inc., USA). HR recording was performed with a chest strap (H10, Polar Electro Oy, Finland) and a wireless receiver (Heart rate Monitor Interface, Sparkfun Electronics, USA) connected to the Simulink model via a serial port. HR measurements were received at a rate of 1 Hz and then downsampled to a sample rate of 0.2 Hz (sample period 5 s) by averaging consecutive batches of five individual samples.
Data preprocessing, model identification, and outcome measures
As noted above, each participant completed two identification tests, thus generating individual data sets (I and II) for model parameter estimation and validation. To prevent overfitting and to eliminate potential orderofpresentation effects, a counterbalanced crossvalidation approach was implemented: for each participant, data set I was used to estimate model parameters and data set II was used as validation data for the estimated models; then, for the same participant, data set II was used for model estimation and data set I for validation. Thus, for the 11 participants, a total of 22 estimation data sets and 22 validation data sets were obtained.
According to the test protocol (Sect. 5.2, Fig. 4a), an evaluation interval from 290 s to 2085 s was used to estimate and validate model parameters. This interval, within one single PRBS period, was selected, such that the number of samples where the input was high (\(v = v_{\mathrm{m}} + 0.25~\)m/s) equalled the number of samples where the input was low (\(v = v_{\mathrm{m}}  0.25~\)m/s). Here, on the evaluation period from 290 s to 2085 s and with a sample period of 5 s, the total number of samples was N = 360, thus giving 180 low samples and 180 high samples.
To remove any potential drifting Phase III dynamic of the HR response, the mean value and any trend were removed (Matlab “detrend” function) prior to estimation and validation; the mean value of the input signal was also removed. An exemplary data set following this preprocessing procedure is provided (Fig. 1), with raw data are shown above (Fig. 4b).
For each estimation data set, two linear timeinvariant transfer functions were employed to model the dynamic response from treadmill speed to HR: a firstorder transfer function (Eq. 3) which combined Phases I and II into a single time constant, and a secondorder transfer function (Eq. 4) with separate time constants for Phases I and II. Hence, for the 11 participants, a total of 22 firstorder models and 22 secondorder models were estimated:
Here, \(k_1\) and \(k_2\) are steadystate gains, and \(\tau _1\), \(\tau _{21}\), and \(\tau _{22}\) are time constants. Model parameters were calculated from the estimation data sets using a leastsquares optimisation procedure (“procest” function from the Matlab System Identification Toolbox; The Mathworks, Inc., USA).
After model estimation, the corresponding validation data sets were used to compute goodnessoffit measures for the resulting first and secondorder models. Two outcome measures were used: the normalised rootmeansquare error [denoted fit, Eq. (5)], and the rootmeansquare error [denoted RMSE, Eq. (6)], as follows:
Here, \(\mathrm{HR}_{\mathrm{sim}}\) is the simulated HR response obtained using the estimated models and the input signal, and HR is the measured HR from the validation data. \(\bar{\mathrm{HR}}\) is the mean value of \({\mathrm HR}\). i is the discrete time index and N is the number of discrete samples considered (as described above, \(N = 360\)). Both of the above outcomes were calculated using the “compare” function from the Matlab System Identification Toolbox.
Statistics
Statistical analysis was performed to test the hypothesis that the goodnessoffit outcomes of secondorder models are better (higher fit and lower RMSE) compared to firstorder models. Prior to analysis, normality of differences between the goodnessoffit outcomes was formally assessed using the Matlab “lilliefors” function (this implements a Kolmogorov–Smirnov test with correction according to the Lillifors method). As it transpired that all differences were not significantly different from a normal distribution, paired onesided t tests were employed for hypothesis testing. Hypothesis testing used a significance threshold of \(5\%\) (\(\alpha = 0.05\)). The Matlab Statistics and Machine Learning Toolbox (The Mathworks, Inc., USA) was employed.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.
Abbreviations
 bpm:

Beats/min
 CI:

Confidence interval
 fit/NRMSE:

Normalised rootmeansquare error
 HR:

Heart rate
 HRV:

Heart rate variability
 MD:

Mean difference
 PRBS:

Pseudorandom binary sequence
 RMSE:

Rootmeansquare error
 SD:

Standard deviation
 k :

Steadystate gain
 \(\tau \) :

Time constant
References
 1.
Hunt KJ, Fankhauser SE. Heart rate control during treadmill exercise using inputsensitivity shaping for disturbance rejection of verylowfrequency heart rate variability. Biomed Signal Process Control. 2016;30:31–42.
 2.
Garber CE, Blissmer B, Deschenes MR, Franklin BA, Lamonte MJ, Lee IM, Nieman DC, Swain DP. American College of Sports Medicine Position Stand. Quantity and quality of exercise for developing and maintaining cardiorespiratory, musculoskeletal, and neuromotor fitness in apparently healthy adults: guidance for prescribing exercise. Med Sci Sports Exerc. 2011;43(7):1334–59.
 3.
Riebe D, Ehrman JK, Liguori G, Magal M, editors. ACSM’s guidelines for exercise testing and prescription. 10th ed. Philadelphia, USA: Wolters Kluwer; 2018.
 4.
Mezzani A, Hamm LF, Jones AM, McBride PE, Moholdt T, Stone JA, Urhausen A, Williams MA. Aerobic exercise intensity assessment and prescription in cardiac rehabilitation. Eur J Prev Cardiol. 2013;20(3):442–67.
 5.
Weston M, Taylor KL, Batterham AM, Hopkins WG. Effects of lowvolume highintensity interval training (HIT) on fitness in adults: a metaanalysis of controlled and noncontrolled trials. Sports Med. 2014;44:1005–17.
 6.
Ramos JS, Dalleck LC, Tjonna AE, Beetham KS, Coombes JS. The impact of highintensity interval training versus moderateintensity continuous training on vascular function: a systematic review and metaanalysis. Sports Med. 2015;45:679–92.
 7.
Hunt KJ, Fankhauser SE, Saengsuwan J. Identification of heart rate dynamics during moderatetovigorous treadmill exercise. BioMed Eng Online. 2015;14:117.
 8.
Wasserman K, Hansen JE, Sue DY, Stringer WW, Sietsema KE, Sun XG, Whipp BJ. Principles of exercise testing and interpretation. 5th ed. Philadelphia, USA: Lippincott, Williams and Wilkins; 2011.
 9.
Whipp BJ, Ward SA, Lamarra N, Davis JA, Wasserman K. Parameters of ventilatory and gas exchange dynamics during exercise. J Appl Physiol. 1982;52(6):1506–13.
 10.
Bearden SE, Moffat RJ. VO2 and heart rate kinetics in cycling: transitions from an elevated baseline. J Appl Physiol. 2001;90(6):2081–7.
 11.
Cheng TM, Savkin AV, Celler BG, Su SW, Wang L. Nonlinear modeling and control of human heart rate response during exercise with various work load intensities. IEEE Trans Biomed Eng. 2008;55(11):2499–508.
 12.
Girard C, Ibeas A, Vilanova R, Esmaeili A. Robust discretetime linear control of heart rate during treadmill exercise. In: Proc. 24th Iran. Conf. Electr. Eng. (ICEE), 2016; pp. 1113–1118.
 13.
Sambeda Sarkar AS. Quasi sliding mode control: an application to heart rate regulation. In: Proc. Int. Conf. Control Power Commun. Comput. Technol. (ICCPCCT), 2018; pp. 299–304.
 14.
Esmaeili A, Ibeas A, Herrera J, Balaguer P, Herrera J, Esmaeili N. Identification and robust control of heart rate during treadmill exercise at large speed ranges. J Control Eng Appl Inform. 2019;21:51–60.
 15.
Du D, Hu Z, Du Y. Model identification and physical exercise control using nonlinear heart rate model and particle filter. In: Proc. 15th Int. Conf. Autom. Sci. Eng. (CASE), 2019; pp. 405–410.
 16.
Baig D, Su H, Cheng TM, Savkin AV, Su SW, Celler BG. Modeling of human heart rate response during walking, cycling and rowing. In: Proc. Ann. Int. Conf. IEEE Eng. Med. Biol. 2010; pp. 2553–2556.
 17.
Argha A, Ye L, Su SW, Nguyen H, Celler BG. Heart rate regulation during cycleergometer exercise using damped parameter estimation method. In: Proc. 38th Ann. Int. Conf. IEEE Eng. Med. Biol. Soc. (EMBC), 2016; pp. 2676–2679.
 18.
Ljung L. System identification: theory for the user. 2nd ed. Upper Saddle River, New Jersey, USA: Prentice Hall; 1998.
 19.
Åström KJ, Wittenmark B. Computer controlled systems: theory and design. 3rd ed. Mineola, New York, USA: Dover Publications; 2011.
Acknowledgements
Lars Brockmann (Institute for Rehabilitation and Performance Technology, Bern University of Applied Sciences) critically reviewed the manuscript for important intellectual content.
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. 320030185351).
Author information
Affiliations
Contributions
KH and HW designed the study. HW did the data acquisition. HW and KH contributed to the analysis and interpretation of the data. HW wrote the manuscript; KH revised it critically for important intellectual content. Both authors read and approved the final manuscript.
Corresponding author
Ethics declarations
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. 201902184). All participants provided written, informed consent.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Wang, H., Hunt, K.J. Identification of heart rate dynamics during treadmill exercise: comparison of first and secondorder models. BioMed Eng OnLine 20, 37 (2021). https://doi.org/10.1186/s12938021008757
Received:
Accepted:
Published:
Keywords
 Heart rate dynamics
 System identification
 Treadmills