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 ([Formula: see text]), O2 uptake ([Formula: see text]), CO2 output ([Formula: see text]) 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 [Formula: see text] and [Formula: see text], whereas for [Formula: see text] and HR the results were similar for both algorithms. Best predictions were found for [Formula: see text] 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, [Formula: see text] and HR responses to incrementally increasing exercise intensities can be sufficiently predicted.

Introduction

The response of cardio-respiratory parameters (CRP), such as oxygen uptake (), ventilation (), CO2 output () 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.

Methods

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). 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 intervals of which were pseudo-randomly assigned to 30 W or 80 W.

Fig. 1

WR protocol, individual measured and predicted response from the calculated impulse response hTD(t) (see text). 1: interval to calculate the mean at 30 W; 2: PRBS interval for hTD(t) calculation; 3: interval to calculate the mean at 80 W; 4: interval to calculate the response to incrementally increansing WR. WR work rate; oxygen uptake. Differences for predictions from hTD(t) and hFD(t) were not visible for this subject

During the exercise test, , and 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 was realized by synchronizing and interpolating all input (WR) and CRP output data (, , , 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 and CCF were derived, using the PRBS samples of the WR as input signal x(t), i. e. WRPRBS(t), and the corresponding CRP data as output y(t), i. e. CRPPRBS(t). As illustrated in Fig. 2, two different approaches were taken to get impulse response functions, and , respectively:

Fig. 2

Obtaining an impulse response function from both, time domain (hTD(t)) and frequency domain (hFD(t)) analysis, respectively. Starting with PRBS changes in WR and the related CRP response, the CCF is calculated and converted into hTD(t)) and hFD(t)

First, in the time domain (TD) the impulse response function was calculated according to Bo et al. (2006):

where and .

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 and . Next, the transfer function was calculated according to the known input–output relation in the frequency domain:

Applying the inverse Fast Fourier transformation to , a second impulse response function was obtained.

For the valid range is limited to due to the composition of the PRBS. Although, is defined for a range , for the following calculations only the range was considered, assuming that under physiological conditions the response is completed after 150 s. To get a comparable basis, was normalized through its sum over the range :

In the last step, the normalized impulse response functions and , must be converted into a relation of WR and the CRP of interest. By convolution of the input with the normalized response function h(t)

and a following linear regression with and the coefficients gain a and offset b are obtained. Therefore, the predictions and as a response to any WR are given for each CRP of interest by:where represents and , respectively, and the regression parameters a and b vary with TD, FD, and each CRP.

Data analysis

The measured data (, , , HR) as well as the predictions and 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 (rPRBS, r80–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 (ResME), the related standard error (ResSE), and the individual correlation coefficients (Res) between Res(x) and measured data were calculated.

Statistical analyses

For the comparisons of rPRBS and r80–75 W resulting from PredTD and PredFD 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 ResME, ResSE and Res with the factors “interval” (PRBS interval and the phase of incremental exercise (80–175 W)) and “parameter” (, , , 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.

Results

Figure 1 shows an individual example for a calculation of in response to the complete exercise protocol.

The CCF data from the PRBS phase and the resulting impulse responses hTD(t), hFD(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), hTD(t) of HR showed the fastest response followed by the responses of and, similarly, and . For hFD(t), different observations can be summarized: hFD(t) of HR started less steep and the maximum for hFD(t) of is much higher compared to hTD(t) of .

Fig. 3

Cross-correlation function (CCF) (a) and the resulting impulse responses hTD(t) (b), hFD(t) (see text for details) (c) for the different cardiorespiratory parameters (means ± SE, N = 13). The grey lines in the background of b/c indicate theoretical responses of first-order systems with time constants τ = 20 s, τ = 30 s and τ = 40 s calculated with these algorithms

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 hTD(t) is zero by definition at t = 0 and t = 150 s, whereas hFD(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 predicted via PredFD. Applying the impulse in the time domain (PredTD), levels at 30 W and 80 W can be predicted with small deviations from measured data (− 3%, 1%).

Table 1

Characteristics for the estimated data and the results from steady-state analysis

Parameter	Method	30 W		80 W		rPRBS	r80–175 W
		mean ± SD	∆	mean ± SD	∆
HR [min−1]	SteadyState	91 ± 11		120 ± 16
	PredTD	95 ± 13	− 4%	115 ± 15	4%	0.64 ± 0.12	0.91 ± 0.09
	PredFD	94 ± 13	− 3%	116 ± 15	4%	0.68 ± 0.11	0.91 ± 0.09
	Significance	a, c		a, b, c		c
V˙_E [L min−1]	SteadyState	21.0 ± 2.5		32.6 ± 2.5
	PredTD	22.2 ± 2.2	− 6%	31.2 ± 2.6	4%	0.45 ± 0.10	0.78 ± 0.14
	PredFD	22.4 ± 1.9	− 7%	31.0 ± 2.7	5%	0.47 ± 0.08	0.78 ± 0.14
	Significance	b		a, b
V˙O_2 [L min−1]	SteadyState	0.71 ± 0.07		1.18 ± 0.08
	PredTD	0.73 ± 0.07	− 3%	1.17 ± 0.09	1%	0.62 ± 0.09	0.82 ± 0.07
	PredFD	0.80 ± 0.07	− 12%	1.08 ± 0.09	8%	0.59 ± 0.09	0.82 ± 0.07
	Significance	b, c		b, c		c	c
V˙CO_2 [L min−1]	SteadyState	0.62 ± 0.08		1.09 ± 0.08
	PredTD	0.66 ± 0.09	− 7%	1.05 ± 0.1	4%	0.57 ± 0.09	0.86 ± 0.06
	PredFD	0.7 ± 0.07	− 13%	1.01 ± 0.08	8%	0.55 ± 0.09	0.86 ± 0.06
	Significance	b, c		b, c		c	c

SteadyState: individual average over 30 s for measured CRP; PredTD: Prediction calculated from hTD(t); PredFD: Prediction calculated from hFD(t); rPRBS: individual correlation coefficient for measured and predicted data of the PRBS interval; r80–175 W: individual correlation coefficient for measured and predicted data of the incremental exercise interval

a: SteadyState vs. PredTD

b: SteadyState vs. PredFD

c: PredTD vs. PredFD

Another criterion to evaluate the algorithm can be derived from the correlations between measured and predicted values (rPRBS, r80–175 W). Significant differences were found for the two prediction algorithms (P < 0.05) for HR (rPRBS only), and (Table 1). The highest correlations for rPRBS were found for HR followed by . showed slightly lower values compared to while 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 r80–175 W were higher than rPRBS in all cases. It can be noted that the ranking between and changed for r80–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 PredTD). The relative residuals (ResME) and standard errors of these (ResSE) (see Eq. 7) differ between HR and , , as well as . 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—PredTD) 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.

Table 2

Analysis of relative residuals (Res) between PredTD and measured data

			HR	V˙_E	V˙O_2	V˙CO_2
ResME	PRBS	x¯	− 0.002	− 0.008	− 0.009	− 0.009
		s	0.001	0.005	0.005	0.005
			a
	80–175 W	x¯	− 0.071	− 0.115	− 0.037	− 0.103
		s	0.049	0.062	0.041	0.055
					a
ResSE	PRBS	x¯	− 0.002	− 0.005	− 0.005	− 0.005
		s	0.001	0.001	0.001	0.001
			a
	80–175 W	x¯	− 0.002	− 0.007	− 0.006	− 0.006
		s	0.001	0.002	0.001	0.002
			a
Resr	PRBS	x¯	0.753	− 0.762	− 0.887	− 0.794
		s	0.100	0.077	0.050	0.072
				a	b
	80–175 W	x¯	− 0.248	− 0.533	− 0.827	− 0.790
		s	0.472	0.219	0.137	0.145
			c		d

Mean (ResME), standard error (ResSE) and correlation (Res) were calculated individually for each CRP in intervals PRBS and 80–175 W. The results were analysed by two-way ANOVA and post-hoc comparisons

aP < 0.05 for all other CRP

b-: p < 0.05

cHR-, HR-: p < 0.05

d-, -: p < 0.05

, mean

s, standard deviation

Fig. 4

Differences ∆ (= PredTD—measured data) for the different CRP related to the rounded measured values for the PRBS phase (black lines) and in the range for 80–175 W (grey lines) (means ± SE, N = 13)

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 (hTD, hFD) 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, hTD and hFD. 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).

showed significant differences between hTD(t) and hFD(t). This difference is the result of the applied normalisation (Eq. 4). While the hTD (t) returns to zero at t = 150 s, this must not be true for hFD(t). For hFD(t) of there was a significant undershoot (see Fig. 3, t > 110 s). This results in a small normalization coefficient and an amplified hFD(t). As stated in the methods section, for hFD(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 hFD(t) in the range 0 < t < 300 s but these analyses yielded non-reasonable and inconsistent results, with higher ResME. Hence, the applied algorithm in the frequency domain must be regarded as less qualified for the prediction (PredFD(t)). The further discussion will be focussed on the PredTD (t) for 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 were found in a similar range with a slightly better predictability of at the 80 W level. rPRBS in Table 1, ResME, and ResSE (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 (r80–175 W) are not surprising due to the increasing CRP values. More interesting are the relative residuals. The significantly lowest mean (ResME) was found for , slightly higher variations were observed for HR and the largest variations were demonstrated for and .

With increasing WR the predictability becomes less accurate, especially for and . 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 hTD(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 hTD(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 , , and 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 PredTD (Table 2) for the 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. and , 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, and . 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 and might be influenced by aerobic capacity and/or aerobic-anaerobic threshold. This might also be true for higher WRs for HR and .

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 and 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, 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.