### 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 (30-min bouts, 3 times per week) and non-smoker, 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 feedback-control 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 vigorous-intensity regime). The mean speed of the treadmill during the final 2 min of the warm up phase was subsequently used as the mid-level 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 fifth-order 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 real-time 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 over-fitting and to eliminate potential order-of-presentation effects, a counterbalanced cross-validation 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 time-invariant transfer functions were employed to model the dynamic response from treadmill speed to HR: a first-order transfer function (Eq. 3) which combined Phases I and II into a single time constant, and a second-order transfer function (Eq. 4) with separate time constants for Phases I and II. Hence, for the 11 participants, a total of 22 first-order models and 22 second-order models were estimated:

$$\begin{aligned} u \mapsto y \negthickspace : \; P_1(s)= & {} \frac{k_1}{\tau _1 s + 1}, \end{aligned}$$

(3)

$$\begin{aligned} u \mapsto y \negthickspace : \; P_2(s)= & {} \frac{k_2}{(\tau _{21} s + 1)(\tau _{22} s + 1)}. \end{aligned}$$

(4)

Here, \(k_1\) and \(k_2\) are steady-state gains, and \(\tau _1\), \(\tau _{21}\), and \(\tau _{22}\) are time constants. Model parameters were calculated from the estimation data sets using a least-squares 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 goodness-of-fit measures for the resulting first- and second-order models. Two outcome measures were used: the normalised root-mean-square error [denoted fit, Eq. (5)], and the root-mean-square error [denoted RMSE, Eq. (6)], as follows:

$$\begin{aligned} \text {fit (NRMSE)} \; [\%]= & {} \left( 1-\sqrt{\frac{\sum _{i=1}^{N}( \mathrm{HR}(i)-\mathrm{HR}_{\mathrm{sim}}(i))^2}{\sum _{i=1}^{N}(\mathrm{HR}(i) - \overline{\mathrm{HR}})^2}} \right) \times 100~\%, \end{aligned}$$

(5)

$$\begin{aligned} \text {RMSE} \; [\mathrm{bpm}]= & {} \sqrt{\frac{1}{N}\sum _{i=1}^{N}(\mathrm{HR}_{\mathrm{sim}}(i) - \mathrm{HR}(i))^2}. \end{aligned}$$

(6)

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 goodness-of-fit outcomes of second-order models are better (higher fit and lower RMSE) compared to first-order models. Prior to analysis, normality of differences between the goodness-of-fit 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 one-sided 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.