Cardiorespiratory kinetics in exercise physiology: estimates and predictions using randomized changes in work rate

Purpose Kinetics of cardiorespiratory parameters (CRP) in response to work rate (WR) changes are evaluated by pseudo-random binary sequences (PRBS testing). In this study, two algorithms were applied to convert responses from PRBS testing into appropriate impulse responses to predict steady states values and responses to incremental increases in exercise intensity. Methods 13 individuals (age: 41 ± 9 years, BMI: 23.8 ± 3.7 kg m−2), completing an exercise test protocol, comprising a section of randomized changes of 30 W and 80 W (PRBS), two phases of constant WR at 30 W and 80 W and incremental WR until subjective fatigue, were included in the analysis. Ventilation (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}_{{\text{E}}}$$\end{document}V˙E), O2 uptake (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}{\text{O}}_{2}$$\end{document}V˙O2), CO2 output (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}{\text{CO}}_{2}$$\end{document}V˙CO2) and heart rate (HR) were monitored. Impulse responses were calculated in the time domain and in the frequency domain from the cross-correlations of WR and the respective CRP. Results The algorithm in the time domain allows better prediction for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}{\text{O}}_{2}$$\end{document}V˙O2 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}{\text{CO}}_{2}$$\end{document}V˙CO2, whereas for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}_{{\text{E}}}$$\end{document}V˙E and HR the results were similar for both algorithms. Best predictions were found for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}{\text{O}}_{2}$$\end{document}V˙O2 and HR with higher (3–4%) 30 W steady states and lower (1–4%) values for 80 W. Tendencies were found in the residuals between predicted and measured data. Conclusion The CRP kinetics, resulting from PRBS testing, are qualified to assess steady states within the applied WR range. Below the ventilatory threshold, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{V}{\text{O}}_{2}$$\end{document}V˙O2 and HR responses to incrementally increasing exercise intensities can be sufficiently predicted.


h TD (t)
Non-normalized Impulse response function resulting from time-domain analysiŝ WR PRBS (t) Approximation of WR in the PRBS interval ϕ xx (jω) Power-spectral density ϕ xx (t) Auto-correlation function ϕ xy (jω) Cross-spectral density

Introduction
The response of cardio-respiratory parameters (CRP), such as oxygen uptake ( V O 2 ), ventilation ( V E ), CO 2 output ( V CO 2 ) and heart rate (HR) to changes in work rate (WR) are measured during cardio-pulmonary exercise tests. From these data, the regulation of aerobic metabolism, ventilation and the cardio-vascular system can be derived (see Poole and Jones 2012 for a comprehensive overview).
In terms of control technology, the CRP responses to exercise can be described as a single-input/single-output system (i. e. WR-CRP). The description of these systems allows to identify the influence of factors such as type of exercise (Koschate et al. 2019a, b), ambient conditions (Drescher et al. 2018) and individual characteristics, e. g. age (Koschate et al. 2016a, b;Ebine et al. 2018;Patti et al. 2021) or fitness (Beltrame et al. 2020;Inglis et al. 2021), on the regulation of the cardiorespiratory system. From the individual kinetics and other characteristics of regulation, interventions, e. g. exercise prescriptions for athletes and patients, can be defined. Furthermore, the combination of exercise tests with information from other non-invasive methods to assess cardiovascular data, e. g. continuous blood pressure measurements, cardiac output and pulse-oximetry, can improve the diagnostic outcome (Oyake et al. 2021). Typically, kinetics are analyzed and described using step responses in combination with data fitting procedures to approximate single mono-exponential functions or sums thereof (Whipp and Rossiter 2005;Keir et al. 2014;Murias et al. 2011). The disadvantage of these methods is twofold: (1) Due to an unfavourable signal-to-noise ratio, several repetitions of the WR steps are required, which extends the test duration significantly and may necessitate more than one test session. This may exclude specific groups of subjects, as older persons, persons with limited endurance, or persons with a tight schedule, from testing.
(2) This kind of modelling requires the assumption of an explicit model with fixed parameters, e. g. a first-order system with time delay (e.g. Ma et al. 2010). The model parameters are determined by iterative least-square criteria and steady states are mandatory. However, if responses to increases (on-step) and decreases in WR (off-steps) are evaluated separately (e.g. Özyener et al. 2001;Fukuoka et al. 2002), potential asymmetries indicate objections to the assumption of dynamic linearity to describe the system response. Changes in step amplitudes may result in differences for the model parameters which is in contradiction to dynamic linearity.
Another approach excites the CRP by WR changes via pseudo-random binary sequences (PRBS) (Hoffmann et al. 2013). This approach requires dynamic linearity between WR and the analyzed CRP, but no explicit model. One PRBS sequence consists of Z intervals that remain constant for a fixed time and are pseudo randomly assigned to changes between two input levels, i.e. WR levels. The resulting sequence with a duration δ Z will be repeated and analysed by time-series analyses in terms of auto-and cross-correlation functions (ACF, CCF). The periodic ACF of the PRBS approximates a periodic impulse function with the periodicity δ Z. In the frequency domain the power spectrum remains flat over a chosen range of frequencies (see Khoo 2018, p. 250ff for a detailed overview).
Z and can be adjusted to concentrate the exciting frequency range to the corresponding frequencies of the CRP dynamics (Seborg et al. 2016, p. 113) and, therefore, the test can be adapted to the given process. For the description of a linear, time-invariant (LTI-) single-input/single-output system the resulting CCF of the input, i.e. WR, and the output, i.e. a CRP, allows a kinetics description in both the time and the frequency domain. The resulting non-parametric description is strongly recommended, when the "…actual model order or time delay is unknown…" and the dynamic behaviour "…cannot be described by standard low-order models" (Seborg et al. 2016, p. 373f).
Both analyses can be summarized by an impulse response function h(t). This function h(t), in combination with any input signal x(t), allows a prediction of the output y(t) by applying the convolution integral for any LTI-system: While h(t) can easily be derived for parametric models, e. g. single mono-exponential functions, the conversion of the resulting CCF from a PRBS test is more complex.
In the following, the application of two different algorithms to convert the CCF of WR-CRP from a PRBS testing into an appropriate impulse response h(t) will be demonstrated. We hypothesize that steady states and the response to incremental WR can be predicted by convolution with the impulse response.
This attempt also reduces the complications caused by the signal-to-noise ratio for the CRP and, by that, reduces the testing time significantly.

Subjects, test protocol and instrumentation
The data for this analysis were taken from the baseline tests of the four missions of campaign four inside the Human Exploration Research Analogue (HERA) facility at NASA Johnson Space Center (see Koschate et al. 2021 for details). (1) Ethical approval was obtained from the Institutional Review Board at the NASA Johnson Space Center (Protocol number: Pro2320) as well as the Ethical Committee of the German Sport University Cologne (Protocol number: 074/2016). Written informed consent was derived from all participants prior to the experiments. 13 subjects (age: 41 ± 9 years, BMI: 23.8 ± 3.7 kg m −2 ) were included in the analysis. The tests were performed on a cycle ergometer in the upright body position, using a complex exercise protocol, which comprised a section with randomized changes between 30 and 80 W (PRBS), two phases of constant WR at 30 W and 80 W and finally an incremental test until subjective fatigue (see Fig. 1). The PRBS consisted of two identical sequences of 300 s. Each sequence was divided into Z = 15 intervals of = 20s which were pseudo-randomly assigned to 30 W or 80 W.
During the exercise test, V O 2 , V CO 2 and V E were measured breath by breath using a metabolic cart (Metalyzer 3B, Cortex Biophysik GmbH, Leipzig, Germany). In addition, beat-to-beat HR was obtained using a wireless ECG belt system (CustoGuard 3, Customed, Ottobrunn, Germany). A constant sampling period t = 1s was realized by synchronizing and interpolating all input (WR) and CRP output data ( V O 2 , V CO 2 , V E , HR). Breath-by-breath data were interpolated stepwise and beat-to-beat data linearly.

Calculation of kinetics response from data in the PRBS interval
The ACF Φ xx (t) and CCF Φ xy (t) were derived, using the PRBS samples of the WR as input signal x(t), i. e. WR PRBS (t), and the corresponding CRP data as output y(t), i. e. CRP PRBS (t). As illustrated in Fig. 2, two different approaches were taken to get impulse response functions,h TD (t) and h FD (t) , respectively: First, in the time domain (TD) the impulse response function h TD (t) was calculated according to Bo et al. (2006): For the second approach, the correlation functions were converted into frequency domain signals, using Fast Fourier transformation. Since the signal power mainly distributes on a few low frequencies for the sequence chosen, the cut-off frequency was determined to be half the PRBS shift frequency (Uhrig 1970, p. 297). Thus, a lowpass filter with passband frequency 25 mHz, which corresponds to the seventh harmonic frequency of this sequence, was applied to the spectral-density functions Φ xx (j ) and Φ xy (j ) . Next, the transfer function H FD (j ) was calculated according to the known input-output relation in the frequency domain:  (t)) and frequency domain (h FD (t)) analysis, respectively. Starting with PRBS changes in WR and the related CRP response, the CCF is calculated and converted into h TD (t)) and h FD (t) Applying the inverse Fast Fourier transformation to H FD (j ) , a second impulse response function h FD (t) was obtained.
For h TD (t) the valid range is limited to 0s ≤ t < 150s due to the composition of the PRBS. Although, h FD (t) is defined for a range 0 ≤ t < 300s , for the following calculations only the range 0s ≤ t < 150s was considered, assuming that under physiological conditions the response is completed after 150 s. To get a comparable basis, h (t) was normalized through its sum over the range 0s ≤ t < 150s: In the last step, the normalized impulse response functions h TD (t) and h FD (t) , must be converted into a relation of WR and the CRP of interest. By convolution of the input WR PRBS (t) with the normalized response function h(t) and a following linear regression with Ŵ R PRBS and CRP PRBS the coefficients gain a and offset b are obtained. Therefore, the predictions Pred TD (t) and Pred FD (t) as a response to any WR are given for each CRP of interest by: where h(t) represents h TD (t) and h FD (t) , respectively, and the regression parameters a and b vary with TD, FD, and each CRP.

Data analysis
The measured data ( V O 2 , V CO 2 , V E , HR) as well as the predictions Pred TD (t) and Pred FD (t) were individually averaged for the last 30 s of the constant WR phases at 30 W and 80 W (refer to Fig. 1). For the PRBS interval and the phase of incremental WR from 80 to 175 W, the Pearson correlation coefficients of measured and predicted data were individually calculated (r PRBS , r 80-175 W ).
Relative residuals Res(t) were individually calculated by: for each second and each CRP. For the PRBS interval and the phase of incremental exercise (80-175 W) individual means (Res ME ), the related standard error (Res SE ), and the individual correlation coefficients (Res r ) between Res(x) and measured data were calculated.

Statistical analyses
For the comparisons of r PRBS and r 80-75 W resulting from Pred TD and Pred FD paired t-tests were applied. Two-way ANOVA for repeated measurements was used to analyse the steady states (factors: WR and method) as well as Res ME , Res SE and Res r with the factors "interval" (PRBS interval and the phase of incremental exercise (80-175 W)) and "parameter" ( V O 2 , V CO 2 , V E , HR). In case of significant effects, post-hoc Bonferroni-tests were applied to compare single means. Level of significance was set to P ≤ 0.05. Figure 1 shows an individual example for a calculation of V O 2 in response to the complete exercise protocol. The CCF data from the PRBS phase and the resulting impulse responses h TD (t), h FD (t) are given in Fig. 3. The time courses for the different CRP differ significantly between the two algorithms (Fig. 3b, c). In contrast to a theoretical first-order model with an instantaneous increase to the maximum and an exponential decrease, the smooth increase (indicated as grey line in Fig. 3b, c), visible for the first seconds (0 < t ≤ 10 s), had to be expected, considering the characteristics of the CCF. The averaging over the individual results with potential time-delays may result in a further smoothing of the signal. In this interval (0 < t ≤ 10 s), h TD (t) of HR showed the fastest response followed by the responses of V O 2 and, similarly, V E and V CO 2 . For h FD (t), different observations can be summarized: h FD (t) of HR started less steep and the maximum for h FD (t) of V O 2 is much higher compared to h TD (t) of V O 2 .

Results
All impulse responses (Fig. 3b, c) calculated from the CCF (Fig. 3a) showed a slow decrease towards t = 150 s. Only for HR a steep decrease can be observed after the maximum until approximately t = 40 s. While h TD (t) is zero by definition at t = 0 and t = 150 s, whereas h FD (t) shows deviations from 0.
The predicted steady-state calculations for 30 W and 80 W were compared with averages from measured data (Table 1). Both algorithms predicted higher 30 W and lower 80 W steady states for all parameters. Remarkably, the highest absolute differences were found for V O 2 predicted via Pred FD . Applying the impulse in the time domain (Pred TD ), V O 2 levels at 30 W and 80 W can be predicted with small deviations from measured data (− 3%, 1%).
Another criterion to evaluate the algorithm can be derived from the correlations between measured and predicted values (r PRBS , r 80-175 W ). Significant differences were found for the two prediction algorithms (P < 0.05) for HR (r PRBS only), V O 2 and V CO 2 ( Table 1). The highest correlations for r PRBS were found for HR followed by V O 2 . V CO 2 showed slightly lower values compared to V O 2 while V E showed the lowest correlations. The wider data range of the CRP data for the increasing WRs compared to the PRBS is the explanation that coefficients r 80-175 W were higher than r PRBS in all cases. It can be noted that the ranking between V O 2 and V CO 2 changed for r 80-175 W .
A more detailed analysis can be derived from the residuals as differences between predicted and measured data ( Table 2 and Fig. 4 for Pred TD ). The relative residuals (Res ME ) and standard errors of these (Res SE ) (see Eq. 7) differ between HR and V E , V O 2 , as well as V CO 2 . For all CRP, the measured data from increasing WR (80-175 W) were significantly underestimated according to the individual means. However, in most cases remarkable correlations between the residual (measured-Pred TD ) were found. The course of residuals (Fig. 4) indicates another aspect for the quality of the predictions: In most cases, systematic decreases of ∆ with increasing values of the respective parameter are visible. Only for HR, a range of stagnation can be identified for lower HR and in the 80-175 W interval.

Discussion
The results of the applied analyses demonstrate the possibility to describe the kinetics of the CRP by an impulse response, and the predictability of CRP steady states as well as responses to increasing WR within certain limits from the PRBS data assessed applying moderate WR intensities. The quality of analysis and prediction must be evaluated separately for each CRP. For all parameters of interest, contradictions to dynamic linear models were identified.

The algorithms
Both analytical approaches (h TD , h FD ) assume dynamic linearity between WR and the respective CRP, since they were calculated from the CCF for the respective CRP. Hence, the mentioned contradictions to the assumptions of dynamic linearity might introduce uncertainties regarding the reliability of the predictions. However, the applied analytical approaches allow to avoid the assumption of a specific model, as used e. g. for exponential fits. However, as demonstrated in Fig. 3, the resulting impulse response may be compared to explicit, e.g. first-order models.
The application of the CCF in conjunction with randomized WR changes of 30 W and 80 W has a disadvantage: The PRBS chosen (15 intervals, 20 s each, total duration 300 s per sequence) has a limited frequency bandwidth. Therefore, rapid changes cannot be analysed. Consequently, in the frequency domain, the analysis must be restricted to the first seven harmonics (period length of 43 s) of the Fast Fourier transformation (Eq. 3) and its inverse (Uhrig 1970, p. 297). This explains the transient increase in the first few seconds of the impulse response for both, h TD and h FD . Even in a theoretical first-order system a mono-exponential decrease, as demonstrated in Fig. 3, cannot be identified. The analysis of simulated responses of first-order systems do not show the typical rapid increase at t = 0 s and the exponential decrease (see also Fig. 3).
VO 2 showed significant differences between h TD (t) and h FD (t). This difference is the result of the applied normalisation (Eq. 4). While the h TD (t) returns to zero at t = 150 s, this must not be true for h FD (t). For h FD (t) of V O 2 there was a significant undershoot (see Fig. 3, t > 110 s). This results in a small normalization coefficient ∑ 149 0h i (t) and an amplified h FD (t). As stated in the methods section, for h FD (t) the impulse response was cut at t > 150 s. The reason for this cut is that for t > 150 s no further physiological origin can be expected as impulse response. However, for some single data sets, calculations were performed with h FD (t) in the range 0 < t < 300 s but these analyses yielded non-reasonable and inconsistent results, with higher Res ME . Hence, the applied algorithm in the frequency domain must be regarded as less qualified for the prediction (Pred FD (t)). The further discussion will be focussed on the Pred TD (t) for V O 2 and the other CRP.

Predictability for the CRP responses
The predictability for the different CRP during moderate exercise, i. e. between 30 and 80 W, can be ranked by means of different criteria and aims. For the steady states, HR and V O 2 were found in a similar range with a slightly better predictability of V O 2 at the 80 W level. r PRBS in Table 1, Res ME , and Res SE (Table 2) indicate the best predictability for HR for the PRBS. The tendency to overestimate lower and to underestimate higher CRP values during the PRBS (Fig. 4) can be assigned to restriction of the analysis to a time frame of 150 s as result from the PRBS applied. The alternative would be an extension of the protocol (Khoo 2018, p. 250). However, the intention was to create a feasible, short test, which can be applied to subjects and patients with restrictions regarding exercise duration and intensity.
The main challenge of this analysis was the prediction of the parameter's behaviour for higher WRs (here 80-175 W). The high correlations (r 80-175 W ) are not surprising due to  1 3 the increasing CRP values. More interesting are the relative residuals. The significantly lowest mean (Res ME ) was found for V O 2 , slightly higher variations were observed for HR and the largest variations were demonstrated for V E and V CO 2 . With increasing WR the predictability becomes less accurate, especially for V E and V CO 2 . This might be related to the increase of anaerobic metabolism and additional ventilatory drives compared to the moderate WR during PRBS. Further research might be focussed on the relation with ventilatory thresholds (see Galán-Rioja et al. 2020 for an overview).

Dynamic linearity
The CRP impulse responses h TD (t) deviate from responses of first-order systems and its time constants as reported in the literature (e.g. Linnarsson 1974;Scheuermann et al. 2002). The shape of the HR h TD (t) seems to be close to a first-order system with a 40 s time constant, even though this would be a much slower response compared with other published data with τ < 30 s (e. g. Linnarsson 1974;Tiedt et al. 1975;DeLorey et al. 2004). For the respiratory parameters, as V O 2 , V CO 2 , and V E show a constant decrease after their maxima which is untypical for first order systems per se. Therefore, these are important arguments to characterize the kinetics of CRP without explicit models.
Although the relative residuals for Pred TD (Table 2) for the V O 2 data indicate impressively low deviations from measured data, this parameter like all others showed significant trends towards higher values (Fig. 4). This is an indication for the need of a non-linear model which will be the subject for further development. The reasons for these non-linearities should be different for each CRP. Commonly, PRBS testing comprises both on-and off-stimuli. The linear modelling implies symmetrical on-and off-responses which might not meet physiological reality. For HR, regulation of blood pressure and vasomotor control are potential physiological explanations for asymmetries. V O 2 and V CO 2 , as pulmonary parameters, might be influenced by non-linearities from venous transport. The portion of anaerobic metabolism leads to additional ventilatory drives with a strong influence on both, V E and V CO 2 . It can be speculated that such non-linearities are caused by the influence from the circulatory system (Eßfeld et al. 1991). Therefore, at least the non-linearities of V E and V CO 2 might be influenced by aerobic capacity and/or aerobic-anaerobic threshold. This might also be true for higher WRs for HR and V O 2 .

Limitations
The PRBS protocol applied in this study should be adapted to the abilities and capacities of the subjects. As demonstrated by Kusenbach et al. (1999) maximal WR in the PRBS protocol may be adapted to specific requirements of the tested subjects. Other modifications were already shown by applications with other modes of exercise, e.g. treadmill (Koschate et al. 2016a, b) or arm cranking exercise (Drescher et al. 2015). With a certain loss of reliability, the PRBS could be shortened, e.g. from 600 to 450 s, which might be a compromise for specific groups.
For V O 2 and V CO 2 analysis, some further improvements might be expected if algorithms are applied to model the alveolar gas exchange more precisely (Koschate et al. 2019a, b).

Conclusions
In summary, the PRBS test, and the resulting kinetics, described by the impulse response, are qualified to describe the kinetics and to assess steady states within the applied WR range. At least below the ventilatory threshold, V O 2 and HR responses to incrementally increasing exercise intensities can be sufficiently predicted from the kinetics derived in the PRBS WR range of 30 W and 80 W. This offers the possibility to assess static as well as kinetics information of subjects with a 10 min moderate exercise test. These data can give important information about the actual state of fitness of the subject which can be used to prescribe physical training or to define its success.
Ethics approval Ethical approval was obtained from the Institutional Review Board at the NASA Johnson Space Center (Protocol number: Pro2320) as well as the Ethical Committee of the German Sport University Cologne (Protocol number: 074/2016) in accordance with the Declaration of Helsinki (including its amendments until 2013).

Consent to participate Written informed consent was available from all participants ahead of the experiments.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.