Determination of heart rate threshold from heart rate kinetics during maximal graded exercise in soccer players

– The main of the present study was to identify the heart rate threshold based on heart rate kinetics during graded maximal exercise in football players. Twenty-six male football players performed a maximal exercise test (Bruce protocol) on a motor-driven treadmill. Oxygen uptake (VO 2 ) and heart rate (HR) were monitored, recorded and resampled at 3.5Hz. The ventilatory threshold (VT), and respiratory compensation (RC), heart rate deflection points (HRDP1 and HRDP2) and heart rate kinetics threshold (HRT) were determined by computerized methods. The heart rate variability (HRV) was assessed in the frequency domain. The HRT averaged 89.9 ± 1.2 % of the VO 2 peak. The HRT showed poor correlations and significant differences compared with HRDP1 (r = 0.46) and VT (r = 0.51), but was not different from, and highly correlated with, HRDP2 (0.98) and RC (0.90). Bland Altman plots showed all athletes into 95% of limits of agreement, and intraclass correlation coefficient showed good agreements between points obtained from HRT compared with HRDP2 (0.96) and RC (0.98). The HRT was highly correlated with HRDP2 and RC, suggesting it could be a marker for cardiorespiratory fatigue.


INTRODUCTION
Since the classical studies by Wasserman`s group, determination of both anaerobic threshold (AT) and peak oxygen uptake (VO 2 peak) have been shown to reflect crucial aspects of cardiorespiratory capacity and have been used to determine the exercise intensity within aerobic training 1,2 . However, in order to assess these markers, a relatively expensive gas exchange apparatus is required. Moreover, during continuous exercise, oxygen kinetics might present different response patterns depending on the nature of metabolic demands. In brief, in exercises performed below AT, oxygen uptake increases initially as a monoexponential function, followed by a plateau phase of approximately 2-3 min (steady state). Conversely, in exercises performed above AT, a triexponential model seems to be more accurate to describe the oxygen kinetics. A slow rise in oxygen uptake (VO 2 ) on the primary component is observed and a steady state phase is less likely to occur 3 .
The fast component of VO 2 kinetics reflects the ability of the cardiorespiratory system to deliver oxygen (O 2 ), as well as the rate at which O 2 is taken up and utilized by exercising skeletal muscle. During moderate exercise, when VO 2 reaches steady-state, the O 2 cost for ATP generation is sustained by aerobic metabolism. During exercises performed with high intensity, the rise of VO 2 (slow component) is probably related to a critical and increased participation of the anaerobic metabolism, which entails progressive lactic acid accumulation 4 .
Previous studies have shown that VO 2 and heart rate (HR) kinetics follow a similar pattern during different exercise modes and intensities 5,6 . Sietsema et al. 6 for instance, have shown that in exercise bouts performed at intensities below AT, HR kinetics present a monoexponential shape. At exercise intensities above AT, the HR would be better described by bi-or triexponential models, with an additional slow rise and no steady-state phase. Thus, a best fit nonlinear least-square approach (migrating from monoexponential to triexponential curves) could be applied to HR kinetics during maximal incremental exercise, and the precise point where the AT occurs 7 .
Therefore, it would be useful to determine whether the transition of HR from monoexponential (with steady-state phase) to either bi-or triexponential fitting (with no steady-state phase) would reflect the metabolic transition associated with AT. This would have practical implications for exercise prescription -if the transition of HR kinetics could be used to determine the AT, it would be possible to establish the intensity of aerobic training without the need of maximum cardiopulmonary exercise testing. Thus, the main goal of the present study was to apply a novel, low-cost and computerized method, based on nonlinear least-square analysis, to identify the HR threshold (HRT) based on HR kinetics in a sample of soccer players. Additionally, the relationship between the HRT determined by this method and AT estimated using ventilatory and heart rate thresholds during maximal incremental exercise has been investigated.

Participants
Twenty-six male football players (age, 26.2 ± 4.8 years; body mass, 73.3 ± 7.8 kg; height, 176.4 ± 7.2 cm) with at least five years of aerobic and resistance training experience volunteered to participate in the study. The experiments were conducted in preseason. All volunteers were considered healthy based on their health history, physical examination, and resting electrocardiogram (ECG). All volunteers answered the Physical Activity Readiness Questionnaire (PAR-Q ) and signed an informed consent according to the Declaration of Helsinki. The experimental procedures were approved by institutional research ethics committee from the university where the study was developed (protocol number 004.2011).

• Ventilatory thresholds
The ventilatory threshold (VT) and respiratory compensation (RC) point were automatically determined by searching for breakpoints of VE/ VO 2 and VE/VCO 2 , respectively. It is known that VO 2 , VCO 2 and VE increase similarly up to VT. However, buffering of lactic acid above VT leads to a disproportional increase in VCO 2 relative to VO 2 with a subsequent increase in VE. Furthermore, the ventilatory equivalent for VCO 2 (VE/VCO 2 ) remains constant or decreases slightly, while the ventilatory equivalent for VO 2 (VE/ VO 2 ) increases. Above the RC point, VE rises at a higher rate than VCO 2 , with a concomitant increase in VE/VCO 2 . Briefly, the breath-by-breath values for VE/ VO 2 and VE/VCO 2 were fitted by the least-squares method to a fifth-degree polynomial using a smooth spline curve. The minima obtained from the first-order derivatives of the fitted polynomials for VE/ VO 2 and VE/VCO 2 were used to calculate the VT and RC values, respectively. Finally, the VO 2peak was computed as the maximal value of VO 2 reached at the end of the test 8 .
The VO 2 and HR were resampled (cubic spline curve) at 3.5 Hz, resulting in equally spaced intervals, and then were filtered by a zero-crossing moving average of length 90 points. Subsequently, data was modeled for each work rate transition for all volunteers, by means of a curvilinear regression obtained from the lowest pooled residual sum of squares.
• Heart rate kinetics, heart rate threshold and heart rate deflection points The work rates were fitted with two models, as follows: model 1, single monoexponential function -HR (t) = BHR + A 1 [1 -e -(t -TD1)/t1and model 2, double exponential function -HR (t) = BHR + A 1 [1 -e -(t -TD1)/ t1 ] + A 2 [1 -e -(t -TD2)/t2 ], where HR (t) is the HR at any time t; BHR is baseline HR (the average HR for the last 30 s of the 3-min warm-up); A 1 and A 2 are HR amplitudes for the fast and slow components, respectively; e represents the exponential fitting; TD1 and TD2 are time delays for the fast and slow components, respectively; and t1 and t2 are time constants for the fast and slow components after their time delays, respectively. The criterion for the HRT determination was set as the last work rate in which steady-state could be reached ( Figure 1). Additionally, the HR time course was modeled by fitting three contiguous linear segments, and two deflection points (HRDPs) were obtained, as previously described 7 . Figure 1. Determination of heart rate threshold (HRT) from a representative subject. Steady-state phase can be observed in the four initial, but not in the last work rate, HRT so being determined at the last work rate (in which steady-state could not be observed).
All heart rate variability (HRV) signal processing was performed using custom algorithms written in Matlab (Mathworks TM , El Segundo, CA, USA). The intervals between heart beats (RR intervals) were extracted at rest and during exercise. For spectral analysis of HRV (frequency-domain), RR intervals were initially resampled at 3.5 Hz, using the spline cubic interpolation method, and the linear trend was removed. Power spectral density was determined by a fast Fourier transform (FFT) based method, the LF (0.045-0.15 Hz) and HF (0.151-1.00 Hz) oscillatory components being reported both as absolute power values (ms 2 ) and as normalized units (n.u.), calculated as follows: HFn.u. = (HF ms 2 )/((LF ms 2 + HF ms 2 ) X 100) and LFn.u.= (LF ms 2 )/((LF ms 2 + HF ms 2 ) X 100).

Procedures
The football players were instructed to eat a light lunch 2 hours before the test. No coffee, tea, or alcohol intake was allowed for 12 h before the experiment and athletes avoided formal and strenuous exercise for at least 48 h before testing. Tests were performed on a motor-driven treadmill (Inbramed TM 10200, Porto Alegre, RS, Brazil) using the Bruce protocol 9 . The football players were allowed sufficient practice time during preliminary testing to become familiar with the treadmill. The tests were preceded by a 3-min warm-up followed by incremental increases in speed and grade every 3 minutes. The following criteria were adopted to consider the test as maximum: a) maximum voluntary exhaustion; b) respiratory exchange ratio R > 1.1; c) maximum HR > 95% of predicted value for age (220 -age); d) VO 2 plateau (VO 2 change ≤ 0.2 l.min -1 between two successive workloads; e) Borg CR10 ≥ 9. Ambient temperature and relative humidity during the tests ranged from 22 to 24°C and 50 to 70%, respectively, and spring water was provided ad libitum.
The VO 2 mask and equipment were put on the subject after they were positioned to perform the exercise, and before the standardized warmup. A face mask (Hans Rudolph TM , Kansas, MO, USA) that covered the mouth and nose of the participant was attached to a bidirectional digital flow valve and fastened by means of a mesh hairnet and Velcro straps. Gas exchange data -VO 2 , carbon dioxide production (VCO 2 ), and VE -were acquired every three complete respiratory cycles, at rest, during, and after treadmill incremental exercise testing via metabolic cart (VO2000, Medical Graphics TM , St. Paul, MN, USA). The ergospirometer was calibrated prior to each individual test according to the manufacturer's instructions. Heart rate was continuously monitored using a beat-to-beat heart rate monitor system (Polar TM RS800, Kempele, Finland).

Statistical analysis
Data normality was tested by D'Agostino & Pearson Omnibus test and all data is presented as mean±sd. Mean values of VO 2 corresponding to HRT were compared to those of VT, HRDP1, RC and HRDP2 using Student's paired t-test, and the correlation between variables was expressed by Pearson's correlation coefficient (r). The agreement between points detected for HRT and HRDP1, VT, HRDP2 and RC were assessed by means of the intraclass correlation coefficient (ICC), and limit of agreement (LOA) was calculated as described by Bland and Altman 10 . The reliability of AT markers were expressed as the differences between VO 2 values obtained by the automatic method and the mean of values detected on the abscissa, which allowed the 95% LOA to be defined as being within 1.96 standard deviations (SD) of the mean difference 10 . Statistical significance for differences between mean values of HRT, VTT, RC and HRDPs were performed using the GraphPad Prism software (version 5). The level of significance was set at p < 0.05. All data were processed using Matlab v. 7.4 R2010a (Mathworks TM , El Segundo, CA, USA).

RESULTS
The cardio-respiratory performance and physical characteristics of athletes are summarized on Table 1.

Relationships between VT, HRDP1, RC and HRDP2 versus HRT
Weak correlations for VT and HRDP1 vs. HRT were found (r = 0.51, p < 0.05 and r = 0.46, p < 0.01, respectively). On the other hand, correlations between RC and HRDP2 vs. HRT proved to be stronger (r = 0.89, p < 0.001 and r = 0.97, p < 0.001, respectively), as well as the correlation between RC and HRDP2 (r = 0.92, p < 0.001). 82.9 ± 6.2 Note. BM: body mass; HR rest: resting heart rate; HR peak: peak of the heart rate; VO 2 peak : peak of the O 2 uptake; HRDP1: first heart rate deflection point; HRDP2: second heart rate deflection point; VT: ventilatory threshold; RC: respiratory compensation point; HRT: heart rate threshold.

Relationships between HRT and frequency-domain of heart rate variability
The power spectral density below HRT showed high predominance of LF compared to HF, but above HRT HF predominates over LF. The results of LF and HF in normalized units (n.u) were obtained at each work rate. Mean values of LF, when compared to HF component, were higher in workloads 1 and 2, but not after the third workload (p < 0.001), when HF increased compared to LF (Figure 3).

DISCUSSION
The present study aimed to apply a novel, low-cost and computerized method, based on nonlinear least-square analysis, to identify the HR threshold (HRT) based on HR kinetics. The major findings of this study are that HRT presented weak correlations with both HRDP1 and VT, but strong correlations with HRDP2 and RC. These results suggest that the present automated computerized method based on nonlinear least-square analysis of the HR time course, may be useful as a marker of cardiorespiratory fitness and fatigue index.
Since the 1960s, several approaches for detection of ventilatory and metabolic thresholds, based either on visual or automated methods, have been developed 1,11 , but some limitations are often discussed. Firstly, gas exchange analysis requires high-cost, sophisticated ergoespirometric equipment and, secondly, the assessment of lactate thresholds (LTs) depends on invasive procedures, consisting of blood collection performed sparsely, which compromises a precise determination of the metabolic transition 6,12 .
After the initial work by Conconi et al. 5 , HR was suggested as a good non-invasive and low-cost method to detect markers associated with fatigue and performance. It should be noted, though, that in the past this approach presented very conflicting results compared to AT, due to low sampling rates at which HR was recorded. Recently, our group developed a computerized regression method to detect performance-associated markers based on HR measured with a relatively high sampling frequency (≥ 1 Hz). HR deflection points (HRDPs) were established -these points did not relate well to AT or VT, but showed high correlation with RC (r = 0.98) 7 . Interestingly, the present results show that the point in which HR kinetics stopped presenting steady-state coincided with the second deflection in HR (HRDP2), which is related to a decrease in HR kinetics amplitude as well as to RC.
Therefore, the issue that arises from these results is: What might be the rationale for VT, HRDP1 and HRT dissociation, and RC, HRDP2 and HRT association? Previous studies have shown the first breakpoint to occur at 44% of VO 2 peak 13,14 , which is within the range of work rates associated to a decrease in the stroke volume (40-50% of VO 2 peak ) and consequently to an increase in afterload. This compensatory HR response seems to occur to stabilize the cardiac output 13,14 . There is no evidence that the decrease in stroke volume has any relationship with HRDP1 or VT, and our results are consistent with the idea that a relationship between these variables likely does not exist. In the present study, the rationale for HRDP 1 after the second work rate (Figure 1) would be a possible association between the increase in HR kinetics amplitude and a higher sympathetic outflow, as denoted by a reduction of HR variability and an increase of low to high frequency ratio (LF/HF ratio) detected by frequency domain analysis.
Some studies suggested that HR increase during incremental exercise would not be linear 7,15 . In this case, it would be very difficult to determine, either by visual or computerized regression methods, the point where AT occurs in graded tests as the Bruce protocol, which concurs with our results. Sietsema and colleagues 6 demonstrated that HR had a slight linear association with VO 2 during graded exercise. Moreover, that HR kinetics presented a monoexponential characteristic at intensities below AT, which became bi-or triexponential at intensities above AT.
During moderate exercise, when steady-state is reached, the O 2 cost for ATP resynthesis is sustained by aerobic metabolism, but during high intensity exercise, the additional rising in VO 2 (slow component) is strongly related to anaerobic metabolism pathways and with blood lactate acummulation 16,17 . Conconi et al. 5,15 suggested that the pattern of HR increase within incremental exercise would not be linear when ATP production relies in a great extent on anaerobic metabolism. In brief, the increase in VO 2 would occur at a lower rate at intensities higher than lower AT, and, as the work rate continues to rise, HR increase seems to be lower compared to VO 2 . Accordingly, our data showed that the point extracted at the transition from mono to bi-or triexponential adjustment during graded exercise, e.g. the HRT, strongly correlated with RC and HRDP2, but not with VT or HRDP1.
With regard to the physiological mechanisms underlying the occurrence of RC, HRDP2 and HRT, they have been classically attributed to an inadequate oxygen delivery and consequent rising in pulmonary ventilation to compensate for the exercise-induced increase in arterial carbon dioxide pressure 1,18 . This higher response of pulmonary ventilation would mostly rely on greater breathing frequency, while tidal volume tends to remain constant 1,19,20 . In parallel, Perlini et al. 19 have shown a concomitant increase in the HF component and breathing frequency during exercise, which has been consistently ratified by further studies 21,22 .
Based on our results, we have presented two possible evidences for the influence of breathing frequency on HRT and HRDP2: first, during rest condition, the HF component oscillates from 0.15 to 0.50 Hz, but, as previously reported 21 , during high intensity exercise, breathing frequency oscillates up to around 1Hz, so it is important to make adjustments in the HF band upper limit (from 0.50 to 1Hz). Second, there was an increase in the normalized HF component with a concomitant decrease in the normalized LF component until work load 3, which could be directly linked to HRDP 2 appearance, since a decrease in the amplitude of HR kinetics was shown above HRDP2.

CONCLUSIONS
The HRT was strongly associated with HRDP2 and RC, but not with HRDP1 and VT. The increase in the HF component of HRV in exercise intensities above HRT might be associated with the mechanical modulatory effect of increased VE on HR. In addition, the present data suggest that this new automated computerized method based on nonlinear least-square analysis of the HR may be a useful tool for exercise training prescription, as well as a marker of cardiorespiratory fitness and fatigue index.

Funding
Silvio Rodrigues Marques Neto is a research fellow at Estácio de Sá University. This study was funded by the Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro -FAPERJ for the "Young Scientist of Our State" fellowship (No. E-26/203.237/2016) granted to Geraldo de Albuquerque Maranhão Neto.

Ethical approval
Ethical approval was obtained from the local Human Research Ethics Committee -Gama Filho University and the protocol (no. 004.2011) was written in accordance with the standards set by the Declaration of Helsinki.